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Neutron  wave  and  pulse  propagation  analysis  is  a  natural  approach 
to  space- dependent  kinetics  in  subprompt  critical  media.  Prior  to  the 
present  work,  analytic  treatments  of  fast  media  have  been  few  and 
limited  in  scope,  in  contrast  to  thermal  wave  and  pulse  propagation, 
which  has  been  studied  rather  thoroughly  and  with  quite  sophisticated 
techniques . 

The  principal  difference  between  analysis  of  fast  and  thermal 
systems  is  treating  the  slowing-down  operator.   A  formal  approach  is 
presented  for  arbitrary  slowing-down  operators;  the  spectrum,  eigen- 
functions,  and  adjoint  eigenfunctions  of  the  slab-geometry  energy- 
dependent  wave  transport  operator  are  obtained,  using  the  singular 
eigenfunction  technique.  Both  multiplying  and  nonmul tip lying  media  are 
treated.  Fission  is  modeled  by  a  one-term  separable  kernel,  although 


Vlll 


the  extension  to  a  multiterm  degenerate  fission  kernel,  representing 
several  fissionable  species,  is  apparent. 

The  fast  neutron  wave  singular  eigenfunction  results  are  compared 
with  other  energy-dependent  transport  work,  particularly  with  previous 
thermal  eigenfunction  analysis,  and  with  work  on  static  fast  neutron 
transport  using  an  energy  transform  approach.  Wave  transport  in  fast 
multiplying  media  and  in  thermal  noncrystalline  media  (modeled  by  a 
separable  thermalization  kernel)  are  rather  similar  in  that  due  to 
energy- regenerative  interaction  processes  a  discrete  asymptotic  separable 
eigenmode  exists  for  moderate  values  of  wave  frequency  and  absorption 
cross  sections.  The  dispersion  laws  obeyed  by  the  eigenvalues  associated 
with  these  modes  are  qualitatively  quite  comparable.  The  fast  nonmulti- 
plying  case  has  no  direct  thermal  analogue  other  than  the  nonphysical 
"absorption  only"  model.   It  is  found  that  the  presence  of  down-scattering 
in  this  case  gives  rise  to  singular  continuum  eigenfunctions  which  are 
not  as  simply  interpreted  as  the  straightforward  streaming  modes 
obtained  for  zero  scattering  cross  section.  Nevertheless,  these  results 
appear  to  be  in  qualitative  agreement  with  other  work  on  energy-dependent 
fast  neuton  transport  theory. 

The  formal  analytic  results  are  developed  in  several  directions  to 
investigate  their  applicability  to  practical  calculations.  A  major  por- 
tion of  this  work  is  devoted  to  obtaining  the  energy-dependent  wave 
transport  representation  of  the  transfer  matrix  method,  which  provides  a 
formalism  for  implementing  calculations  concerning  wave  and  pulse  propa- 
gation through  finite  regions  such  as  adjacent  slabs  of  different  compo- 
sition.  It  is  found  that,  for  isotropic  scattering,  the  basic  operators 


IX 


of  the  transfer  matrix  formalism  can  be  constructed  from  eigenfunctions 
of  the  wave  transport  operator.   This  result  is  general  and  is  equally 
applicable  to  fast  and  thermal  analysis. 

Finally,  the  dispersion  law  expression  for  a  fast  multiplying 
medium  is  employed  to  develop  an  algorithm  for  computing  the  discrete 
eigenfunctions  and  associated  dispersion  law  for  separable  and  degener- 
ate fission  kernels.   A  specific  application  of  this  method  then  is 
made  to  the  case  of  isotropic  elastic  and  inelastic  scattering  from  any 
number  of  nuclides  and  levels,  with  arbitrary  lethargy  dependence  of 
cross  sections.  Elastic  scattering  is  modeled  by  a  free  gas  kernel,  and 
inelastic  scattering  by  a  constant  energy  loss  per  interaction  per 
level.   Unlike  techniques  requiring  inversion  of  matrices,  computation 
time  increases  approximately  linearly  with  increases  in  lethargy  steps, 
making  quite  detailed  computations  feasible.   Illustrative  computations 
are  carried  out  using  constant  cross  sections  and  a  single  elastic 
scattering  species. 


CHAPTER  I 


INTRODUCTION 


Purpose 

As  the  demand  for  energy  increases  and  supplies  of  economically 
recoverable  fuels  diminish,  fast  breeder  reactors  will  supply  an 
increasing  proportion  of  base- load  generated  power.  To  operate  most 
economically  these  reactors  will  tend  to  be  as  large  as  is  technologi- 
cally feasible.  Accordingly,  as  the  size  of  fast  reactor  cores 
increases,  it  will  become  increasingly  important  to  understand  spatially 
dependent  kinetic  effects  in  fast  systems. 

A  particularly  straightforward  method  of  investigating 
propagation  of  neutronic  disturbances  in  fast  reactor  materials  is  to 
place  a  pulsed  or  oscillating  source  of  neutrons  at  the  face  of  an 
experimental  assembly,  and  then  to  observe  the  propagation  of  the  neu- 
tron "signal"  through  the  assembly.   In  this  way  one  can  study  spatially 
dependent  flux  oscillations  such  as  might  be  expected  to  result  from 
flow- induced  vibrations  of  core  components,  void  formation  and  collapse, 
and  other  such  phenomena. 

Neutron  wave  and  pulse  propagation  experiments  have  been  performed 
in  many  different  thermal  media,  both  multiplying  and  nonraul tip lying. 
The  theoretical  basis  of  analysis  of  thermal  wave  propagation  in 


nonmul tip lying  media  has  attained  a  considerable  degree  of  sophistication, 
and  fairly  accurate  numerical  prediction  of  some  experimental  results  is 
possible.   This  is  in  contrast  to  the  situation  in  the  fast  neutron  wave 
regime;  few  experiments  have  been  performed,  and  analytic  investigations 
have  been  hampered  by  difficulties  which  do  not  arise  in  treatment  of 
thermal  systems.  The  purpose  of  this  dissertation  is  to  present  a 
particular  framework  of  approach  within  which  these  difficulties  may  be 
addressed,  extending  techniques  which  have  been  applied  primarily  to 
thermal  analysis. 

Three  general  objectives  will  be  pursued: 

(i)  to  develop  the  spectral  representation  of  the  energy- 
dependent  fast  wave  Boltzmann  operator  as  far  as  possible 
in  sufficiently  general  form  so  that  its  potential  for  use 
with  realistic  cross  section  data  can  be  evaluated; 
(ii)  to  extend  a  formalism  which  treats  neutron  transport  in 

finite  and  discontinuous  media  so  that  the  above  results  may 
be  applied  to  wave  transport  in  experimentally  realistic 
geometries  and  through  successive  regions;  and 
(iii)  to  illustrate  applications  of  the  analysis  by  computing  the 
fundamental  eigenfunction  and  dispersion  law  for  wave 
propagation  in  fast  multiplying  media,  using  a  modelled 
kernel  in  the  Boltzmann  operator. 

Early  Neutron  Wave  Investigations 

In  1948  Weinberg  and  Schweinler  published  the  first  description  in 
the  open  literature  of  the  generation  and  analysis  of  neutron  waves  [1]. 


Using  one-speed  diffusion  theory  they  were  able  to  show  that  a  localized 
oscillation  in  neutron  absorption  within  a  reactor  would  produce  a 
perturbation  in  the  neutron  population  which  would  propagate  in  wave- 
like fashion.  The  first  experiments  with  neutron  waves  were  reported 
in  1955  by  Raievski  and  Horowitz  [2],  using  a  mechanically  modulated 
exterior  source  to  generate  waves  in  DO  and  graphite.  Uhrig  [3]  then 
applied  this  technique  to  measurements  in  subcritical  assemblies.  Both 
experimental  and  theoretical  aspects  of  neutron  wave  propagation 
subsequently  received  considerable  attention  and  refinement,  particularly 
by  Perez  [4]  and  his  associates  at  the  University  of  Florida,  although 
experiments  and  most  analytic  efforts  were  restricted  to  thermal  systems. 

As  investigation  of  the  theoretical  basis  of  neutron  wave 
experiments  proceeded,  it  was  realized  that  from  an  analytic  standpoint 
experiments  involving  spatially  propagating  pulses  were  equivalent  to 
neutron  wave  experiments,  since  any  physically  realizable  pulse  could  be 
time-Fourier  analyzed  to  give  its  frequency  components  [5] .  Also,  it 
became  clear  that  neutron  wave  propagation  was  related  to  other  linear 
static  and  kinetic  experimental  techniques,  in  particular  the  classical 
exponential  experiment,  which  is  the  zero-frequency  limit  of  the  wave 
experiment,  and  the  pulse  die-away  experiment  [S-8] .   (The  die-away 
experiment  monitors  the  time-rate  of  decay  of  a  neutron  population 
which  has  been  introduced  into  a  finite  assembly  by  a  pulsed 
external  source.   For  technical  reasons  this  type  of  experiment 
was  easier  to  perform  than  wave  experiments,  and  enjoyed  a  more  rapid 
initial  development  [4].)  As  a  result,  these  methods  experienced 


considerable  parallel  theoretical  treatment  [8-10] .  Early  work  in  this 
area  is  reviewed  extensively  by  Uhrig  [7]  and  Perez  and  Uhrig  [4] . 

Fast  Neutron  Wave  Investigations 

Neutron  wave  and  pulse  propagation  has  received  proportionately 
very  little  attention  in  the  fast  energy  regime.  The  only  experiments 
described  in  the  literature,  performed  by  Napolitano  et  al.  [11,12]  and 
Paiano  et  al.  [13]  at  the  University  of  Florida,  have  been  in  nonmulti- 
plying  media;  no  experiments  in  multiplying  media  have  been  reported. 
The  technical  difficulty  of  such  experiments  probably  has  contributed 
both  to  the  lack  of  experimental  data  and  the  scarcity  of  methods  to 
predict  and  correlate  results.  Theoretical  analysis  also  has  been 
retarded  by  the  fact  that  even  tractable  energy-dependent  analytic 
models  of  fast  media  do  not  have  convenient  mathematical  properties, 
and  consequently  most  of  the  elegant  techniques  which  have  been  applied 
to  thermal  neutron  transport  cannot  be  extended  readily  to  this 
problem  [14].  Notable  exceptions  to  the  general  absence  of  numerical 
techniques  and  results  are  the  multi-group,  multiplying  medium  calcula- 
tions of  Travelli  [15,16],  and  the  calculations  of  Booth  et  al.  [17], 
using  the  multi-group  discrete  ordinates  method  of  Dodds  et  al.  [18]  to 
interpret  Napolitano' s  experimental  results. 

Energy- Dependent  Transport  Formulation 
of  Wave  and  Pulse  Propagation 

Before  discussing  the  various  theoretical  results  which  are 

directly  or  indirectly  applicable  to  fast  neutron  wave  and  pulse 


propagation,  it  will  be  helpful  to  approach  the  general  neutron  wave 
problem  from  the  point  of  view  of  the  energy-dependent  transport  method 
which  will  be  used  in  this  dissertation.  We  begin  with  the  classic 
time-dependent  Boltzmann  equation  for  the  neutron  flux  [19,20,21], 
which  we  will  write 

^^  4>(r,E,^,t)  +  ti  '   V(j)(r,E,J^,t)  +  Z^CE)4)(r,E,J^,t) 


d^ 

4lT    0 


dE"  K(E^f^'  ^  E,Q)<piT.E\^\t)   =  S(r,E,fi,t)     (1.1) 


In  standard  notation,  (})(r,E,J^,t)  is  the  neutron  directional  flux,  v  is 
the  neutron  speed,  u   is  the  unit  vector  in  the  direction  of  neutron 
travel,  E  is  the  neutron  energy,  r  is  the  spatial  coordinate,  ^^(E)  is 
the  total  cross  section,  and  S(r,E,U,t)  is  a  flux- independent  source 
term.  The  interaction  kernel  K(E'',S^'  ^  E,ti)   contains  all  neutron  inter- 
action processes  which  give  rise  to  secondary  neutrons  with  altered 
energy  and  direction.  All  cross  sections  (and  thus  K)  are  taken  to  be 
region-wise  independent  of  r,  and  the  medium  described  by  the  equation 
is  assumed  to  be  isotropic  in  the  sense  that  interaction  properties  do 
not  depend  on  the  initial  direction  of  neutron  travel.  Cross  sections 
also  are  assumed  to  be  constant  with  respect  to  time  and  independent  of 
the  flux;  neutron  wave  and  pulse  propagation  experiments  meet  these  two 
criteria  quite  well.  Propagation  of  disturbances  in  sub-prompt  critical 
reactors  also  should  be  adequately  described  by  this  linear  kinetic 


model.   (Nonlinear  space-dependent  kinetics  are  of  interest  primarily  in 
the  context  of  excursion  situations;  such  problems,  while  important,  are 
difficult  to  analyze,  and  thus  far  have  been  approached  by  use  of 
specialized  and  involved  computational  techniques  [22,23,24].) 

Any  of  a  number  of  classic  analytic  approaches  can  be  taken  to  the 
solution  of  Eq.  (1.1);  here  we  will  treat  it  as  an  eigenvalue  problem. 
This  will  enable  us  to  extend  our  results  to  finite  medium  and  multi- 
region  problems;  the  transfer  matrix  formalism,  which  will  be  discussed 
in  Chapter  III,  requires  solution  of  a  similar  eigenvalue  equation,  and 
we  will  be  able  to  relate  its  solution  to  those  of  the  wave  Boltzmann 
equation  for  isotropic  scattering.  Furthermore,  we  can  make  use  of 
spectral  analysis  which  already  has  been  done  on  the  thermal  neutron 
version  of  Eq.  (1.1). 

Plane  Symmetry  and  the  Eigenvalue  Equation 

The  mathematical  development  of  transport  theory  has  reached  its 
greatest  sophistication  for  the  case  of  plane  symmetry,  and  this  is  triie 
also  for  the  particular  subject  of  neutron  wave  and  pulse  propagation. 
Since  this  geometry  also  is  appropriate  for  the  description  of  classi- 
cal wave  and  pulse  propagation  experiments,  we  will  turn  our  attention 
to  the  specific  case  of  plane  neutron  waves.  The  infinite  medium  plane 
wave  eigenfunctions  which  will  be  obtained  may  then  be  used  in  develop- 
ing the  corresponding  transfer  matrix  formalism,  which  can  be  employed 
to  study  the  propagation  of  waves  and  pulses  through  finite  slabs  and 
successive  slabs  of  dissimilar  materials. 


Following  customary  arguments,  we  stipulate  that  all  sources  or 
initial  fluxes  must  be  rotational ly  symmetric  about  the  x  axis  and  do 
not  depend  on  the  transverse  Cartesian  coordinates  y  and  z.   Orienting 
the  X  axis  along  the  direction  of  wave  propagation,  we  have  for  the 
transport  operator  of  Eq.  (1.1) 


^  '   h  =  Mj^  ;  ]s  E  n   '   X  (1.2) 


where  x  is  the  unit  vector  in  the  x  direction,  and  y  is  the  cosine  of 
the  angle  between  the  path  of  neutron  travel  and  the  x  axis.   With  the 
restrictions  of  Eq.  (1.2)  the  homogeneous  Boltzmann  equation  may  be 
written 


:^-^  (})(x,E,y,t)  +  u  ^  4)(x,E,y,t)  +  Z^(E)  (t)(x,E,y,t) 


1     r°° 
dy'   dE'  K(E',y'  -^  E,y)(})(x,E',y',t)  =  0.        (1.3) 

; 

-1      0 


We  notice  that  time  and  space  operators  appear  only  in  the  first 
two  terms,  respectively,  while  the  integral  operator  acts  on  E  and  y. 
Consequently,  x  and  t  variables  may  be  separated.  With  appropriate 
choices  for  the  separation  constants  <))  may  be  expressed  as  a  damped 
plane  wave 


(t)^(x,E,y,t)  =  F(E,y;K)  e^*^^  e"'^'' 


„,„  ,   .   -ax  i(wt-?x) 
=  F(E,y;ic)  e    e  ^     ^ 


(1.4) 


where  u  is  the  wave  or  Fourier- component  frequency,  and  k  is  a  complex 
constant,  K  =  a  +  iC,  the  complex  inverse  relaxation  length.  Thus  a  is 
the  inverse  relaxation  length  of  the  wave,  while  C  is  its  wave  number. 
The  frequency  O)  will  be  regarded  as  a  parameter  of  the  equation,  and  we 
will  treat  K  as  the  eigenvalue  to  be  determined.   Introducing  Eq.  (1.4) 
into  Eq.  (1.3) 


—  -  y<  +  Z^(E) 


F(E,y;K)   - 


dy^ 


dE'  K(E^y-'  ^  E,y)F(E^y^K)    =   0, 

(1.5) 


or  defining 


a(E,a.)    EZ^(E)    .  ^^^^ 


(1.6) 


Eq.    (1.5)   has   the  form 


fl  f°° 


(a(E,a))   -  yK)F(E,y;K) 


dy         dE     K(E    ,y     ->  E,y)F(E    ,y   ;k)    =   0 


-1  o 


(1.7) 


which  will  be  referred  to  as  the  wave  Boltzmann  equation,  or  WBE, 
throughout  the  rest  of  this  work.   This  is  the  most  general  statement  of 
the  Boltzmann  equation  in  wave  eigenvalue  form  with  plane  symmetry.  As 
an  eigenvalue  problem  it  should  more  properly  be  written 


f  1    r°° 


dy^ 


dE' 


^^^  6(E-E')6(y-yO  -  i  K(E^y'  ->  E,y) 


F(E',y^K) 


=  KF(E,y;K).  (1.8) 

The  pole  at  y  =  0  causes  no  difficulties  which  we  will  need  to 
consider  [25] ;  this  value  of  y  corresponds  to  a  direction  of  neutron 
travel  perpendicular  to  the  direction  of  wave  propagation. 

Interaction  Operators  for  the  Fast  and 
Thermal  Neutron  Regimes 

The  two  general  types  of  neutron  interaction  which  are  of 

importance  in  wave  and  pulse  propagation  and  which  enter  into  the 

kernel  of  Eq.  (1.7)  are  scattering  and  fission.  Adopting  for  a  moment  a 

theoretician's  perspective  on  reality,  we  may  define  a  fast  neutron 

experiment  as  one  in  which  the  scattering  kernel  has  a  Volterra  form  in 

energy.   In  a  similar  vein  a  thermal  problem  may  be  distinguished  by  the 

presence  of  a  Fredholm  scattering  kernel.  These  observations  stem  from 

the  fact  that  in  the  fast  neutron  regime  one  is  concerned  with  neutron 

energies  from  the  eV  range  to  about  10  MeV,  the  upper  end  of  the  fission 

spectrum;  hence  only  downscattering  in  energy  is  important.   By  contrast. 
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in  the  thermal  regime  neutrons  are  in  or  near  thermal  equilibrium  with 
their  surroundings  so  that  upscattering  in  energy  occurs  as  well; 
energies  of  interest  range  essentially  over  the  thermal  Maxwellian 
spectrum. 

Appropriate  interaction  kernel  models  reflect  these  properties. 
We  may  write  the  thermal  scattering  contribution  to  the  interaction 
kernel  as 


1 


dy' 


dE'  Z^(E^y"  ^  E,y) 


-1 


while  the  fast  scattering  operator  has  the  form 


■1    f°° 

dy'  I  dE'  Z^(E^y'  ^  E,y) 

1     E 


extending  the  notation  of  Eq.s  (1.3)  to  (1.7). 

In  multiplying  media  the  interaction  kernel  contains  a  contribution 
due  to  fission  in  addition  to  scattering.  Difficulties  associated  with 
treating  the  slowing  down  of  fission  neutrons  [26]  have  precluded  trans- 
port analysis  of  wave  propagation  in  thermal  multiplying  media,  although 
other  models  such  as  age-diffusion  have  been  employed  [27] .  No  such 
problem  arises  in  the  fast  multiplying  wave  problem,  since  the  energy 
range  of  the  fission  spectrum  is  essentially  the  energy  range  of 
interest. 
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For  the  fission  contribution  to  the  interaction  operator  we  will 
use  the  customary  isotropic  separable  kernel 

When  only  one  fissionable  isotope  is  present  this  is  a  satisfactory 
model.  For  two  or  more  fissionable  species,  one  can  either  construct  an 
equivalent  separable  kernel  with  averaged  x  and  vZ-  or  employ  a  degene- 
rate kernel 


j  =  l 


Only  the  separable  kernel  will  be  treated  in  detail  here  but  the 
formalism  of  Chapter  II  can  be  extended  in  a  straightforward  way  to 
degenerate  kernels  for  multiple  fissioning  species. 

To  avoid  the  appearance  of  unwelcome  factors  of  1/2  in  connection 
with  the  isotropic  fission  spectrum,  we  will  make  the  following 
notational  distinction.  Define  the  isotropic  x(E)  so  that 


X(E)dE  =  1.  (1.9) 


Define  X  so  that 


X  =  x(E,y)  =  1/2  xCE);  (1.10) 


12 


then 


.00 


dy 


dE  X  =  1. 


(1.11) 


•1     0 


The  fission  interaction  operator  employing  the  separable  kernel  model 
and  Eq.  (1.10)  thus  becomes 


1    /'°° 

du' 

J 
■1  0 


dE'  vZ^(E')' 


which  is  the  form  which  will  be  used  throughout  this  work. 

Using  the  above  forms  for  contributions  to  the  interaction  kernel 
the  fast  homogeneous  WBE  of  Eq.  (1.7)  may  be  written 


Ca-viK)F(E,ja;K)  = 


/•I 

dy' 
}  i 

1     E 


dE'  E^(E',u'  ^  E,y)F(E',y  ;k)    (1.12) 


for  nonmultiplying  media,   and 


(a-yK)F(E,y;K) 


fl         r°° 
dy' 


dE'  Z    (E',y' 


E,y)F(E',y';K) 


■1  E 


•1  f°° 


+  X 


dy'       dE' 


vZ_(E')F(E',y';K) 


(1.13) 


■1  o 


for  multiplying  media. 
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It  should  be  noted  that  we  have  taken  into  account  only  prompt 
neutrons  in  Eq  (1.13),  and  hence  v  here  is  the  number  of  prompt  neutrons 
per  fission.   Delayed  neutrons  will  contribute  only  at  wave  periods 
greater  than  the  shortest  delayed  neutron  precursor  lifetime,  an  effect 
which  has  been  investigated  numerically  by  Travelli  [16].  Also,  it  will 
be  assumed  that  the  medium  under  consideration  is  subprompt  critical. 

The  thermal  nonmultiplying  WBE,  which  will  be  discussed  as  a  point 
of  departure  for  our  work  on  Eq.  (1.13),  is 


(a-u<)F(E,y;K)  = 
-1 


'1     f°° 

du'  dE'  Z  (E',m'  ^  E,m)F(E',m';k).    (1.14) 


Spectrum  and  Eigenfunctions  of  the 
Thermal  Transport  Operator 

The  fast  neutron  wave  energy- dependent  transport  eigenvalue 
problem  can  best  be  introduced  by  discussing  work  which  has  been  done  on 
the  analogous  thermal  problem,  Eq.  (1.14).  This  approach  will  be  taken 
because  transport  treatment  of  the  fast  problem  is  necessary  to  obtain 
qualitatively  correct  spectral  descriptions  for  passive  media.  Approxi- 
mations such  as  diffusion  theory  can  yield  an  estimate  of  the  least 
attenuated  mode  of  propagation  in  fast  multiplying  media,  where  such  a 
fundamental  mode  exists,  but  can  provide  little  other  information  rele- 
vant to  the  properties  one  should  expect  of  the  exact  transport  treat- 
ment. 

The  context  of  the  work  to  be  presented  here  is  the  "singular 
eigenf unction  method,"  which  received  its  major  impetus  from  a  paper  by 
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Case  [28],  and  thus  is  frequently  known  as  "Case's  method."  As  an 
introduction  to  the  literature  on  the  singular  eigenfunction  method  in 
transport  theory,  including  the  wave  problem,  the  review  of  McCormick 
and  Kuscer  [29]  is  highly  recommended,  as  it  is  both  recent  and 
extensive. 

Travelli  [30]  was  the  first  investigator  to  arrive  at  an 
essentially  correct  description  of  the  spectrum  of  the  energy-dependent 
wave  problem,  based  on  a  multigroup  transport  formulation.  We  turn  now 
to  the  energy- dependent  analysis  of  the  thermal  wave  eigenvalue  problem, 
corresponding  to  Eq.  (1.14),  performed  independently  by  Kaper  et  al.  [31] 
and  Duderstadt  [32,33];  the  former  study  employs  an  isotropic  one-term 
degenerate  thermalization  kernel,  while  Duderstadt  discusses  more  gene- 
ral types  of  scattering  interaction  models  as  well.  Their  results  are 
summarized  in  this  section.  Eq.  (1.8)  may  be  written  in  abbreviated 
form  as 


AF  E  (A  +  A  )F  =  <F  (1-15) 


where  the  streaming  operator  A  is 


rl 

mQO 

dy' 

dE' 

J 

A,  E    du'   dE'  -^ +   ""^ 


1  j  "^  J   I  y   yv(E)j 


6(E"-E)6(y"-y)-     (1.16) 
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and  the  interaction  operator  A  becomes 


^2  = 


fl     j-oo 


dE'  -  I    (E^y' 
y  s 


E,y) 


(1.17) 


-1 


using  the  scattering  interaction  kernel  of  Eq.  (1.14). 

The  basic  method  for  obtaining  the  spectrum  and  eigenfunctions  of 
this  equation  is  a  generalization  of  the  work  of  Bednarz  and  Mika  [34] 
on  the  static  Boltzmann  operator,  which  in  turn  extended  the  classic 
monoenergetic  singular  eigenfunction  technique  [25]  to  a  continuous 
energy  representation.  We  begin  by  defining  the  domain  C  in  the  spectral 
K-plane,  which  is  the  continuous  spectrum  of  the  streaming  operator  A  : 


^^(E) 


1  0) 

yv(E) 


K  =  0,  y  e  [-1,1],  E  e  [o,«>) 


(1.18) 


or  in  the  notation  of  Eq.  (1.14),  those  values  of  <  for  which  a  -  yK 
vanishes.  For  any  nonzero  frequency  co,  a  is  complex,  so  that  C  will 
occupy  an  area  in  the  K-plane.   It  is  instructive  to  consider  both  the 
rectangular  and  polar  forms  of  <  6  C;  Eq.  (1.18)  implies  that 


^JE) 


Re(K) 

=   a   : 

y 

Im(K) 

=  C  = 

0} 

yv(E) 

r(<) 

y 

*-             V   ■' 

e(K) 

=   tai 

-1     w 

1/2 


(1.19) 
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where  r  and  6  are  the  usual  radial  and  azimuthal  polar  coordinates.   In 
general  C  will  consist  of  two  symmetric  portions  in  the  first  and  third 
quadrants,  due  to  u  6  (0,1]  and  y  6  (0,-1]  respectively.  This  is  represented 
schematically  in  Figure  1.1.   Since  frequency  w  is  a  positive  quantity, 
C  does  not  extend  into  the  second  and  fourth  quadrants.   Figure  1.2 
shows  the  first  quadrant  of  the  <-plane  in  more  detail.  The  domain  C 
is  bounded  by  the  line  y  =  ±  1,  a  =  E  ;  from  the  rectangular  form  of 
Eq,  (1.19)  it  is  apparent  that  the  real  part  of  this  boundary  line 
assumes  every  value  of  Z  as  E  (and  thus  v)  varies  from  0  to  °°.  The 
polar  form  of  Eq.  (1.19)  shows  that  as  |y|  varies  from  1  to  0,  values  of 
K  corresponding  to  a  fixed  E  generate  a  line  of  constant  9  which  begins 
at  the  boundary  of  C  and  extends  to  infinity.   As  the  parameter  o)  is 
increased  or  decreased  the  domain  C  expands  or  contracts  in  the  imaginary 
K  direction. 

We  note  that  if  vZ  varies  monotonically  with  E,  each  point  of  the 
domain  C  will  correspond  to  a  unique  E,y  pair,  E  and  y  ;  for  the  ther- 
mal  analysis  presented  here,  this  is  assumed  to  be  the  case.   Two 
important  results  then  follow.   First,  Eq.  (1.19)  defines  a  one-to-one 
mapping  of  E,y  onto  the  spectral  plane.  Second,  the  spectrum  of  A  is 
not  degenerate.   The  consequences  of  these  results  will  be  discussed 
later. 

The  discrete  and  residual  spectra  [35]  of  A  are  empty  [33] .   The 
singular  continuum  eigenfunctions  of  A  ,  satisfying  the  equation 


(A^-K)F(E,y;K)  =  0  (1.20) 
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esIm(K) 


a=  Re(K) 


Figure  1.1    The  continuum  domain  C  in  the  spectral 
K-plane. 
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Figure  1.2    Structure  of  the  continuum. 
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are 


FCE,y;K)  -  6(E-E^)6(u-u^)  =  6(a-y<)  .  (1.21) 


This  gives  a  corresponding  eigenmode,  using  Eq.  (1.4)  (for  A  =  A.), 

<t>    (x,E,]A,t)   =   e'^""  e^*^'^^'^''^6(E-E^)6(M-y  )         (1.22) 

which  clearly  represents  neutrons  of  energy  E  streaming  in  the 
direction  y  .  Since  C  is  in  one-to-one  correspondence  with  all  possible 
E,y  pairs,  each  point  in  the  continuum  corresponds  to  a  unique  neutron 
speed  and  direction  of  travel.  Referring  to  Eq.  (1.19)  we  see  that 
modes  with  y  =  ±  1  have  relaxation  length  1/E  (E  )  equal  to  the  neutron 
relaxation  length;  that  is,  modes  corresponding  to  k  on  the  boundary  of 
C  represent  neutrons  streaming  along  the  x  axis.  Other  modes  are  more 
attenuated,  as  the  direction  of  neutron  travel  becomes  more  oblique  to 
the  direction  of  wave  propagation. 

The  spectrum  of  the  streaming  operator  and  its  eigenfunctions  are 
qualitatively  the  same  for  both  thermal  and  fast  regimes,  the  only  dif- 
ferences being  the  values  of  E  which  are  applicable,  and  the  detailed 
structure  of  E^  as  a  function  of  energy;  Eqs .  (1.19)  -  (1.22)  apply  in 
either  case.   It  is  the  interaction  operator  A  ,  containing  the  descrip- 
tion of  the  scattering  and  multiplication  processes,  which  gives  rise  to 
the  qualitative  differences  between  fast  and  thermal  WBE  eigenfunctions. 

It  seems  likely  for  "reasonable"  mathematical  models  of  thermal 
scattering  that  the  spectrum  of  A  =  A  +  A  always  contains  the  spectrum 
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of  A^.  This  has  been  substantiated  for  A-  having  the  form  of  a  separable 
kernel  with  isotropic  scattering  [31,33].   This  model,  which  was  pro- 
posed by  Corngold  et  al.  [36],  has  been  used  quite  extensively  in 
analytic  transport  studies,  since  it  represents  fairly  well  the  quali- 
tative features  of  thermal  scattering  interactions  [37] .  Duderstadt  [33] 
has  investigated  more  general  operators  A  ,  and  while  the  spectral 
analysis  for  less  restricted  models  is  somewhat  more  tentative,  it  does 
appear  to  indicate  that  the  spectrum  of  the  streaming  operator  is  in 
general  contained  in  the  spectrum  of  the  wave  Boltzmann  operator  A.  We 
will  see  this  in  a  more  formal  way  from  the  technique  used  to  construct 
the  continuum  eigenf unctions. 

To  illustrate  this  method  we  obtain  the  eigenfunctions  for  the 
thermal  WBE  with  the  separable  isotropic  thermal izat ion  kernel 


-1     O 


dE'      ^gCE')-  , 


.00 


du 


J     J 
1    o 


dE  M(E)  =  1.  (1.23) 


(To  be  consistent  with  our  later  treatment  of  the  fast  WBE  we  will  not 
perform  the  usual  symmetrization  of  this  kernel,  since  for  the  fast  case 
A  is  not  symmetrizable.   The  main  result  of  interest  which  arises  from 
symmetry  of  A  is  that  the  eigenfunctions  are  mutually  orthogonal,  and 
one  avoids  the  adjoint  problem;  this  and  other  considerations  will  not 
be  of  direct  concern  here.  Also  note  that  M(E)  is  not  the  Maxwellian 
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distribution.  To  satisfy  detailed  balance  M(E)  =  M'(E)Z  (E) ,  where 
M^(E)  is  proportional  to  the  Maxwellian,  subject  to  the  above  normaliza- 
tion constraint.) 

Using  this  form  for  A  ,  the  WBE  corresponding  to  Eq.  (1.14) 
becomes 


(a-yK)F(E,y;K)  =  M(E) 


dy' 


dE'  Z  (E')F(E',y';K).     (1.24) 


Notice  that  this  expression  is  exactly  equivalent  to  the  fast  multiplying 
WBE  in  the  form  of  Eq.  (1.13)  when  scattering  is  ignored  in  that  equa- 
tion. First  we  investigate  the  point  spectriim.  We  see  that  k  will  be 
an  eigenvalue  when  the  homogeneous  equation  (1.24)  has  a  solution  for 
that  value  of  k.      Let  us  suppose  that  k  0  C  so  that  (a  -  y<)  ^   0;  then 
dividing  by  this  factor. 


F(E,y;K) 


M(E) 
a-yK 


dy- 


dE'  Z  (E')F(E',y'*;<) 
s 


-1 


Defining  the  scalar  product 


(1.25) 


((|)(E,y),  i)iE.]x))   = 


.00 


dy 


dE  (})(E,y)iJ;(E,y) 


(1.26) 
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Eq.  (1.25)  may  be  written  more  compactly 


F  = 


M 


a-yK  ^  s 


(^.,F) 


(1.27) 


Taking  the  scalar  product  of  this  equation  with  E  and  eliminating  the 
scalar  factor  (I  ,F)  we  find  that  the  condition  for  Eq.  (1.26)  to  have  a 
solution  is 


Z   -^ 
s'   a- UK 


(1.28) 


Defining  the  dispersion  function 


A(K,w) 


s'    a-iJK 


(1.29) 


Eq.  (1.28)  is  simply  the  condition  that  this  dispersion  function  vanish. 
Eq.  (1.28),  which  is  referred  to  as  the  dispersion  law,  determines  in 
the  present  problem  the  regular  eigenvalues  k  of  the  WBE  as  a  parametric 
function  of  frequency.   Indeed,  values  of  k  which  satisfy  the  dispersion 
law  for  a  given  frequency  o)  have  been  shown  [31,33]  to  comprise  the 
point  spectrum  of  A  with  A^  defined  by  Eq.  (1.23);  for  these  eigen- 
values K.  the  corresponding  eigenfunctions  are  given  by  Eq.  (1.27): 


F(E,y;<.)  =  A(K.)  ^^^ 
*•  '  '  j'^    ^  ]     a-yK 


(1.30) 
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where  A(k  )  is  an  arbitrary  constant.   (Note  that  for  small  co  and  Z 

■'  a 

this  approaches  a  Maxwellian  distribution  in  energy.) 

When  K  e  C,  the  term  (a  -  y<)  is  zero  for  a  particular  E  =  E  and 

K 

^^  =  \-     As  we  have  mentioned,  C  is  contained  in  the  continuous  spectrum. 
In  the  present  case  the  continuous  spectrum  of  A  is  identically  the 
domain  C,  and  the  continuum  eigenfunctions  are  [31,33,34] 

'^^'^'^^  =W^\^^^  \   ^E^  ^,^^')n^',V^;K)   .   A«)6Ca-pK) 

-1     o 

(1.31) 

K  e  c 

using  the  notation  of  Eq.  (1-21);  A(k)  is  an  arbitrary  constant.  We  see 
that  (a  -  yK)"  has  a  pole  at  the  "eigen-energy"  E^  and  "eigen-angle" 
y^;  integrals  over  E,y  involving  this  term  will  exist  in  the  ordinary 
sense,  provided  that  its  coefficients  in  the  integrand  are  well-behaved 
at  the  pole.  Hence  we  may  eliminate  the  scalar  (Z  ,F)  in  Eq.  (1.31)  in 
favor  of  the  constant  A  (k)  by  taking  the  scalar  product  of  the  equation 
with  Z^(E)  and  solving  for  (Z  ,F).  We  then  find 


F(E,y;K)  =  A(k) 


M(E)  ^s ^^<^ 

a-y< 


J —  +  (S(a  -  yK) 

(1.32) 


K  e  c 


so  that  A(K)  is  in  fact  a  normalization  constant;  A  is  the  dispersion 
function  defined  in  Eq.  (1.29) 
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Eq.  (1.31)  may  be  obtained  directly  from  Eq.  (1.24)  by  a  heuristic 
argument  [38].  Since  for  any  variable  x  the  function  x6(x)  is  identi- 
cally zero,  apparently  X(a  -  yic)6(a  -  y<)  E  0  may  be  added  to  the  right- 
hand  side  of  Eq.  (1.24).   Division  by  (a  -  ijk)  gives  Eq.  (1.31)  when 
K  e  C.  Evidently,  then,  the  domain  C  always  will  be  in  the  spectrum  of 
A,  since  it  is  contributed  by  the  streaming  operator,  regardless  of  the 
form  of  A  . 

The  continuum  eigenfunction,  Eq.  (1.32),  is  composed  of  two 
singular  terms,  one  being  the  pure  streaming  mode  of  Eq.  (1.21),  and  the 
other  having  distributed  E  and  y  dependence,  with  the  same  formal  struc- 
ture as  the  discrete  eigenfunction,  Eq.  (1.30),  except  that  it  has  a 
pole  singularity  since  k  6  C.  The  scalar  coefficient  of  the  latter 
term,  E  (E  )/A,  represents  the  relative  excitation  of  the  distributed 
portion  of  the  mode  by  the  streaiming  portion  (this  can  be  seen  more 
clearly  by  comparing  with  the  analogous  fast  continuum  eigenfunction, 
which  will  be  developed  in  Chapter  II).  Hence  the  continuum  eigenfunc- 
tion may  be  interpreted  as  being  due  to  direct  streaming  neutrons  having 
energy  and  direction  E  and  y  ,  and  an  associated  scattered  distribution 
which  is  excited  by  the  streaming  neutrons;  the  scattered  distribution 
is  peaked  at  E  and  y  due  to  the  pole  of  the  transport  coefficient 
(a  -  yK)   ,  but  contains  all  other  E,y  values  as  well.  Note,  however, 
that  the  entire  mode  has  the  phase  velocity  v   H  ^  (cf.  Eq.  (1.4))  of 
the  uncoil ided  wave.  This  interpretation  of  the  thermal  continuim 
eigenf unctions  has  not  been  given  in  previous  treatments,  as  symmetriza- 
tion  of  the  kernel  tends  to  obscure  the  physics  involved. 
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Kaper  et  al .  [31]  have  investigated  the  dispersion  law  A  =  0  for 

the  separable  kernel  model.  Their  findings  may  be  summarized  as 

follows.  When  03  =  0,  one  has  the  classic  exponential  experiment  [21]; 

there  is  exactly  one  pair  of  real  eigenvalues  k.  =  ±  k  ,  provided  that 

J      P 
the  absorption  of  the  medium  is  not  too  strong  (of  course  the  precise 

condition  will  depend  on  the  energy  dependence  of  the  cross  sections) . 
Otherwise,  the  point  spectrum  is  empty  and  will  remain  empty  for  all  co. 
As  the  parameter  w  is  increased  from  zero,  the  pair  of  eigenvalues  will 
move  symmetrically  into  the  first  and  third  quadrants  of  the  complex 
K-plane.  Evidently  for  sufficiently  large  w  there  will  be  a  limiting 
frequency  u^  beyond  which  the  discrete  spectrum  is  empty;  this  value  of 
frequency  appears  to  occur  when  K  meets  the  boundary  of  the  continuum 
C.  This  situation  is  represented  schematically,  for  the  first  quadrant, 
in  Figure  1.3.   (We  noted  that  in  general  the  boundary  of  C  is 
frequency- dependent;  here  for  simplicity  it  is  shown  for  Z  constant,  in 
which  case  the  boundary  remains  a  line  perpendicular  to  the  real  axis.) 

While  for  a  time  it  was  conjectured  [32,33]  that  zeroes  of  the 
dispersion  function  might  exist  within  the  continuum  as  "embedded  eigen- 
values" in  a  continuation  of  the  dispersion  law  for  w  >  oj  ,  it  now 

c 

appears  [31,39]  that  this  is  not  the  case,  although  the  dispersion 
function  apparently  does  vanish  at  points  within  the  continuum  [39]; 
referring  to  Eq.  (1.32),  this  corresponds  to  points  at  which  the  delta- 
function  contribution  vanishes.  This  subject  will  not  be  pursued  here; 
the  interested  reader  is  referred  to  Kaper  et  al.  [31],  Kline  and 
Kuscer  [39],  and  for  an  extensive  discussion  from  a  different  point  of 
view,  to  the  work  of  Doming  and  Thurber  [40]  and  Doming  [41]. 


26 


le 


/ 


/ 


/ 
/ 


/ 


I   /'^=0 


Figure  1.3    Schematic  dispersion  law  for  a  discrete 
eigenvalue. 
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From  the  results  of  calculations  based  on  the  separable  kernel  or 
comparable  models  [31,32]  it  is  possible  to  draw  some  conclusions  about 
the  physical  interpretation  of  the  frequency- dependent  behavior  of  the 
dispersion  law.   Near  zero  frequency  the  fundamental  neutron  wave  mode 
(if  there  is  one)  corresponding  to  ±  <  is  less  attenuated  than  the 
streaming-associated  continuum  modes.  As  cj  increases,  the  wavelength  of 
the  fundamental  mode  decreases  and  it  becomes  more  attenuated.  This 
occurs  because  it  becomes  increasingly  difficult  for  scattered  neutrons 
to  remain  in  phase  with  the  wave;  from  Eq.  (1.24)  we  see  that  as  < 
approaches  C  the  energy  and  angular  distributions  of  the  wave  become 
increasingly  peaked  for  y  =  ±  1.  Ultimately  the  fundamental  mode 
becomes  nearly  as  attenuated  as  forward- streaming  neutrons,  and  evi- 
dently the  distributed  E,y  term  of  the  continuum  eigenfunctions  then 
assumes  the  role  formerly  held  by  the  fundamental  mode  as  the  frequency 
increases  beyond  w  . 

One  additional  remark  should  be  made.  For  o)  =  0  the  spectrum  of 
the  static  Boltzmann  equation  lies  entirely  on  the  real  axis,  and  in 
general  it  appears  that  it  is  the  limit  of  the  spectrum  of  the  WBE  as 
u)  approaches  zero.  But  obviously  for  k  real  the  identification  of  E,u 
pairs  with  points  of  C  no  longer  can  be  made.   Indeed  it  may  be  improper 
to  regard  the  static  Boltzmann  equation  as  the  zero- frequency  limit  of 
the  WBE,  No  such  problem  arises  in  connection  with  the  discrete  spectrum, 
in  the  sense  that  in  the  limit  00=  0,  Eq.  (1.28)  gives  the  correct 
eigenvalues  for  the  static  case.  This  evidently  is  true  of  the  disper- 
sion law  in  general,  and  in  that  sense  we  speak  loosely  of  the  exponen- 
tial experiment  being  the  zero- frequency  limit  of  the  wave  experiment. 
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In  this  work  we  will  be  concerned  only  with  co  ^  0  except  in  the 
calculations  of  Chapter  IV,  which  involve  only  the  discrete  eigenvalues 
and  eigenfunctions. 

Completeness  of  the  Thermal  Eigenfunctions 

In  order  to  make  use  of  a  set  of  eigenfunctions  such  as  those 
obtained  in  the  previous  section,  it  is  necessary  to  show  that  arbitrary 
functions  (suitably  restricted)  can  be  expanded  using  these  functions 
as  a  basis,  and  it  is  further  necessary  to  evaluate  the  expansion 
coefficients.  First,  then,  one  must  prove  that  the  set  of  eigenfunc- 
tions is  complete,  or  at  least  establish  completeness  within  the  context 
of  the  problem  one  is  to  consider.  Then  either  the  eigenfunctions  must 
be  shown  to  be  orthogonal  and  normalized  to  unity  scalar  product,  so 
that  orthogonality  properties  may  be  used  to  obtain  expansion  coeffi- 
cients in  the  usual  way,  or  some  other  procedure  must  be  followed. 
Normalization  of  continuum  eigenfunctions  is  somewhat  less  than  straight- 
forward because,  as  may  be  seen  from  the  form  of  these  eigenfunctions  in 
Eq.  (1.32),  it  involves  products  of  delta  functions  of  complex  variables. 
The  alternative  procedure  is  to  find  the  continuum  expansion  coeffi- 
cients G(k)  of  an  arbitrary  function  4'(E,u)  directly  from  the  expression 
for  the  expansion,  which  is  a  singular  integral  equation: 


r(E,y)  = 


G(K)  F(E,ii;K)dK  (1.33) 

K 
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where  'F'  is  the  portion  of  ^   contributed  by  the  continuum  eigenfunctions 
(the  discrete  mode  contribution  is  found  by  the  usual  application  of 
orthogonality).  F(E,y;<),  which  is  now  the  kernel  of  the  integral 
operator,  is  known  from  Eq.  (1.32)  or  a  similar  evaluation  of  the  con- 
tinuum eigenfunction  based  on  another  model. 

In  implementing  either  method  the  theory  of  generalized  analytic 
functions  [42]  has  been  the  principle  tool  applied  to  date.  Eq.  (1.33) 
also  has  been  used  to  prove  completeness  of  the  WBE  eigenfunctions, 
since  if  it  can  be  shown  that  an  arbitrary  V(E,\i)   has  a  representation 
in  this  form,  the  set  of  eigenfunctions  F(E,y;K)  must  be  complete. 
This  approach  has  been  taken  by  Kaper  et  al .  [31]  and  Duderstadt  [32] 
to  show  completeness  for  the  eigenfunctions  of  the  separable  thermali- 
zation  model  of  the  previous  section;  their  treatments  were  based  on 
extension  of  the  generalized  analytic  function  technique  as  applied  by 
Cercignani  [43]  to  problems  in  the  kinetic  theory  of  gases.  The  details 
of  this  analysis  are  lengthy  and  will  not  be  repeated  here. 

We  will  make  reference  to  two  types  of  completeness  and 
orthogonality.  We  note  that  values  of  k  in  the  first  quadrant  correspond 
to  plane  waves  propagating  in  the  positive  x  direction,  and  similarly 
the  third  quadrant  represents  waves  propagating  in  the  negative  x  direc- 
tion.  In  general,  e.g.,  within  a  slab  of  finite  thickness,  a  wave  will 
be  made  up  of  components  traveling  in  both  directions;  to  represent  an 
arbitrary  wave  (or  pulse  frequency  component)  '!'(£, y,x,u))  in  WBE  eigen- 
functions, one  must  use  all  the  eigenfunctions,  corresponding  to  the 
whole  spectrum  of  the  wave  Boltzmann  operator.  Completeness  of  the 
first  type,  in  the  sense  that  a  unique  representation  of  this  sort  can 
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be  made,  is  termed  full-range  completeness.  The  corresponding  full- 
range  orthogonality  is  simply  orthogonality  under  the  scalar  product  of 
Eq.  (1.26). 

In  applying  eigenfunction  techniques  to  boundary  value  problems, 
one  frequently  wishes  to  represent  an  incoming  source  flux  S(E,y,x  ,03) 
or  to  specify  flux  continuity  for  waves  moving  from  one  region  into 
another  across  an  interface  at  a  boundary  point  x  .   In  this  case  the 
boundary  condition  will  be  specified  for  either  ye  [-1,0)  ory  G  (0,1]  and 
will  involve  eigenfunctions  for  only  one  direction  of  wave  propagation. 
Completeness  in  this  sense,  termed  half -range  completeness,  requires 
that  a  function  defined  over  y  6  (0,1]  or  y  6  [-1,0)  can  be  represented 
liniquely  by  WBE  eigenfunctions  corresponding  to  the  eigenvalues  in  only 
the  first  or  third  quadrant  of  the  spectral  plane,  respectively.  Half- 
range  orthogonality  is  orthogonality  under  integration  the  half  range  of 

y. 

Both  full-range  and  half-range  completeness  requirements  will  be 
seen  to  arise  in  Chapter  III  in  connection  with  a  formalized  treatment 
of  the  slab  geometry  boundary  value  problem.  We  should  note  that  at 
present  half-range  completeness  can  be  proved  only  for  quite  restricted 
kernel  models,  although  full-range  completeness  can  be  shown  for  more 
general  kernels  [29,33].  Our  main  interest  in  the  completeness  proper- 
ties of  the  eigenfunctions  of  Eqs .  (1.24)  and  (1.26)  is  that  they  are 
indeed  complete.   We  will  use  the  same  formal  procedure  to  find  the 
eigenfunctions  of  the  fast  WBE,  and  will  obtain  qualitatively  similar 
results.  Thus  we  may  have  considerable  confidence,  in  lieu  of  proof, 
that  the  fast  eigenfunctions  are  complete  as  well. 
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There  are  two  difficulties  which  will  prevent  us  from  extending 
the  generalized  analytic  function  technique  directly  to  the  fast  regime. 
First,  one  must  deal  with  the  slowing-down  operator.  Second,  the  one-to- 
one  equivalence  between  values  of  E,y  and  points  of  C  does  not  hold  for 
realistic  fast  cross  sections  (e.g.,  at  resonances),  and  we  will  be 
reluctant  to  consider  more  restrictive  cross-section  models  (i.e., 
monotonic  vE  );  this  equivalence  plays  a  central  part  in  the  generalized 
analytic  function  method  as  it  has  been  developed  to  date.  Whether 
these  two  problems  are  insurmountable  is  a  matter  for  further  investiga- 
tion; however,  it  seems  unlikely,  in  view  of  the  results  established  in 
the  thermal  case,  that  the  fast  eigenf unctions  would  not  be  complete 
for  "reasonable"  cross- section  models.   (An  example  of  an  "unreasonable" 
model  is  a  strictly  1/v-dependent  cross  section 


or  one  which  has  this  behavior  over  some  energy  range.  When  this  occurs 
the  portion  of  C  corresponding  to  this  energy  range  collapses  onto  a 
line.  This  case  is  discussed  for  thermal  waves  in  polycrystalline 
material  by  Duderstadt  [33]  and  by  Yamagishi  [44];  it  is  necessary  to 
deal  separately  with  the  eigenfunctions  on  the  line  continuum  which 
results  from  this  cross  section.) 

For  an  introduction  to  other  literature  on  completeness  of  singular 
eigenfunctions  see  the  review  of  McCormick  and  Kuscer.   It  is  interest- 
ing in  this  connection  to  read  the  comments  of  Burniston  et  al .  [45],  and 


32 


the  recent  remarks  of  Zweifel  [46]  regarding  the  degree  to  which  the 
rigorous  mathematical  basis  for  the  singular  eigenfunction  transport 
analysis  has  been  established. 

Other  Related  Problems  and  Literature 

In  the  foregoing  discussion  we  have  seen  that  for  the  thermalization 
model  employed  there  the  spatially  dominant  wave  mode  is  due  to  the 
regular  eigenvalue  k  ,  which  is  determined  by  the  zeroes  of  the  disper- 
sion function  A(K,aj).   Further,  we  see  from  Eq.  (1.32)  that  the  zeroes 
or  near- zeroes  of  A  also  will  play  a  large  part  in  determining  the 
character  of  continuum  modes,  since  in  regions  where  A  is  small  the 
scattering  portion  of  the  mode  will  dominate  the  streaming  term.  A 
corresponding  dispersion  function  appears  to  arise  in  general  in  the 
treatment  of  regenerative  media  (i.e.  those  in  which  neutron  interac- 
tions can  result  in  either  a  gain  or  loss  in  energy,  and  hence  the 
interaction  kernel  has  a  Fredholm  form) .  Doming  and  Thurber  [40] ,  for 
example,  find  that  in  an  alternative  formulation  of  the  wave  problem  and 
in  an  initial  value  problem  the  nature  of  solutions  is  similarly 
influenced  by  the  behavior  of  a  dispersion  function.   In  addition, 
dispersion  laws  are  known  to  arise  in  nontransport  approximations  to 
dynamic  eigenvalue  problems.   For  example,  when  the  multi-group  diffu- 
sion approximation  is  used  to  obtain  a  matrix  expression  analogous  to 
Eq.  (1.8),  its  determinant  is  the  dispersion  function,  and  the  disper- 
sion law  is  simply  the  requirement  that  the  determinant  vanish;  the 
solutions  associated  with  values  of  k  which  satisfy  the  dispersion  law 
are  then  the  desired  eigenmodes.   Indeed  the  multi-group  diffusion 
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approach  has  been  used  rather  extensively  to  compute  dispersion  laws  for 
moderators,  and  when  sufficiently  accurate  scattering  matrices  are 
employed,  agreement  of  diffusion  theory  methods  with  experiment  at  low 
frequencies  can  be  quite  good  [47]. 

General  discussions  of  wave  and  pulse  propagation  in  the  context  of 
its  relationship  to  other  dynamic  problems,  properties  of  the  various 
dispersion  laws,  and  analytic  methods  which  have  been  applied  to  these 
problems  will  be  found  in  Bell  and  Glasstone  [21]  and  Hetrick  [48].  An 
excellent  review  of  the  literature  in  this  area  as  of  1967  has  been 
given  by  Kuscer  [49],  although  it  is  of  course  somewhat  dated.  As  an 
alternative  exact  approach  to  transport  problems,  the  Wiener-Hopf  tech- 
nique is  finding  increasing  favor  and  must  be  viewed  as  a  potential 
method  for  analysis  of  the  wave  problem;  Williams  [50]  recently  has 
published  an  expository  review  of  the  method.  Also,  the  singular  eigen- 
function  method  review  of  McCormick  and  Ku^Ser  [29]  should  be  mentioned 
again  in  connection  with  the  subject  of  transport  treatments  of  various 
static  and  dynamic  problems. 

Finally,  with  respect  to  the  subject  of  Chapter  III,  note  should  be 
taken  of  existing  work  treating  neutron  waves  in  geometry  which  is 
finite  or  has  discontinuities  along  the  direction  of  wave  propagation. 
Interface  effects  first  were  investigated  experimentally  by  Denning, 
Booth  and  Perez  [51].  This  same  problem  was  the  subject  of  both  numeri- 
cal and  analytic  investigation  by  Baldonado  and  Erdmann  [52,53];  their 
work  is  of  particular  interest  because  one- speed  and  energy-dependent 
diffusion  and  transport  results  are  given.  Mockel  [54]  has  presented 
both  transfer  matrix  and  invariant  imbedding  transport  formulations  for 
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wave  transmission  and  reflection  from  a  slab  imbedded  in  an  infinite 
medium  of  different  composition.  Also  to  be  noted  is  the  treatment  of 
Larson  and  McCormick  [55]  of  transport  in  a  slab,  in  the  static  case, 
using  a  degenerate  scattering  kernel.  Recently  much  attention  has  been 
given  by  Japanese  and  Indian  groups  to  the  problem  of  thermal  neutron 
wave  propagation  in  assemblies  of  polycrystalline  moderating  materials 
(e.g.  graphite  and  beryllium)  having  finite  transverse  and  longitudinal 
dimensions;  see  for  example  Nishina  and  Akcasu  [56],  Kumar  et  al.  [57], 
and  Yamagishi  [44] .  The  latter  is  of  particular  interest  because  it 
demonstrates,  in  a  transport  treatment,  the  presence  of  intermodal 
interference. 


CHAPTER  II 


SPECTRUM  AND  EIGENFUNCTIONS  OF  THE 
FAST  WAVE  TRANSPORT  OPERATOR 


Introduction 

In  this  chapter  the  singular  eigenfunction  formalism,  presented  in 
Chapter  I,  will  be  extended  to  the  fast  WBE  expressions,  Eqs.  (1.12)  and 
(1.13).  Both  the  forward  and  adjoint  eigenfunctions  will  be  obtained 
for  general  forms  of  the  nonmul tip lying,  or  "slowing-down,"  and  multiply- 
ing cases.  The  structure  of  these  solutions  will  be  discussed,  and  some 
of  the  implications  of  using  realistic  cross-section  and  scattering  ker- 
nel models  will  be  explored. 

Adjoint  eigenfunctions  will  be  investigated  for  two  reasons. 
First,  they  will  be  necessary  for  the  treatment  of  the  transport  formu- 
lation of  the  transfer  matrix  in  Chapter  III.  Second,  as  has  been  men- 
tioned, analytic  evaluation  of  expansion  coefficients  cannot  be  performed 
using  generalized  analytic  function  techniques  which  have  been  applied 
to  thermal  problems.  For  the  same  reason,  we  will  not  obtain  normaliza- 
tion constants  analytically.  However,  biorthogonality  of  eigenfunction 
sets  will  be  shown  in  the  classic  way. 

Adjoint  Eigenfunction  Equations 

The  appropriate  scalar  product  under  which  to  define  adjoint 
operators  is  given  by  Eq.  (1.26).  We  consider  the  general  wave  eigenvalue 
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equation  in  the  form  of  Eq.  (1.8),  which  we  may  write  as 


^FCE,u;k)  -1^ 


rl 

1    o 


r°° 


dE'  K(E',u  ->  E,y)F(E^y  ;k)  =  KF(E,y;K) 


(2.1) 


The  adjoint  eigenfunctions  will  be  denoted  by  F"  (K;E,y);  the  adjoint 
eigenvalue  equation  corresponding  to  the  forward  equation,  Eq.  (2.1),  is 


^F"^(K;E,y) 


dy' 
-1     0 


dE"  ^  K(E,y  -^  E^y^)F"^(K;E^y")  =  KF"^(K;E,y) 


(2.2) 


where  we  notice  the  factor  —  is  now  within  the  integral.  However,  if 
we  define 


f''"  (K;E,y)  E  i  F"^(K;E,y) 


(2.3) 


Eq.  (2.2),  becomes,  upon  substitution  and  rearranging. 


(a  -  y<)F"''  (<;E,y)  - 


dy'   dE'  K(E,y 


E^y')F■^(<;E^y^) 


=  0     (2.4) 


■1     o 


which  is  the  form  which  would  have  been  obtained  as  the  adjoint  of  the 
homogeneous  wave  Boltzmann  equation,  Eq.  (1.7).   It  will  be  more  conven- 
ient to  deal  with  Eq.  (2.4)  since  it  differs  from  the  forward  WBE, 
Eq.  (1,7),  only  in  the  kernel  of  the  interaction  operator  and  hence  we 
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will  be  able  to  apply  the  same  techniques  to  the  solution  of  both 
forward  and  adjoint  equations. 

It  should  be  pointed  out  that  k  is  used  for  the  eigenvalue  in 
Eq.  (2.2),  with  the  implication  that  the  spectra  for  forward  and  adjoint 
equations  are  identical.  That  this  is  true  for  "well-behaved"  operators 
in  the  models  we  are  considering  will  be  apparent  from  the  singular 
eigenfunction  formalism,  although  of  course  each  case  must  be  explained 
individually.  Nicolaenko  [14]  has  exhibited  an  inelastic  scattering 
operator  for  which  the  adjoint  spectrum  contains  additional  contribu- 
tions due  to  a  singularity  of  the  kernel  at  zero  energy;  he  uses  the 
singular  kernel  in  defining  an  energy  transform  for  reduction  of  the 
static  transport  slowing-down  equation  (for  the  model  he  considers)  to 
monoenergetic  form.  However  it  is  shown  in  Appendix  B  that  singularity 
of  inelastic  scattering  kernels  is  not  an  inherent  attribute  of  fast 
neutron  transport  problems.  Thus  for  the  forward  and  adjoint  problems 
the  spectra  and  eigenfunctions  can  be  regarded  tentatively  as  being  in 
correspondence,  subject  to  verification  for  specific  interaction  models. 

Biorthogonality  of  Eigenfunctions 

Biorthogonality  of  WBE  eigenfunctions  corresponding  to  different 

eigenvalues  can  be  shown  by  the  usual  argument.  Writing  Eq.  (2.4)  for 

K  and  Eq.  (1.7)  for  k,  we  take  scalar  products  of  the  two  equations 

f 
with  F(K)  and  F  (k')  respectively,  and  subtract  to  find 


(<  -  k') 


yF'''(K'),F(<) 


=  0,  (2.5) 
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noting  that  we  are  using  a  real-type  scalar  product,  Eq.  (1.26).  We 
conclude  that  biorthogonality  holds  for  F  and  F  under  a  u-weighted 
scalar  product,  while  in  view  of  Eq.  (2.3)  this  is  equivalent  to 
biorthogonality  of  F   and  F  with  unit  weighting. 


Nonmultip lying  Media:  Spectrum  of  the 
Slowing-Down  Transport  Operator 

In  fast  nonmultip lying  media  the  WBE  is  given  by  Eq.  (1.12).   The 

corresponding  adjoint  WBE,  Eq.  (2.4),  is  found  to  be 


(a  -  u<)F'''(K;E,y)  = 


fl 


dy' 


dE'  E^(E,y  ^  E^u")F■^(K;E^y')    (2.6) 


where  the  different  energy  limits  for  the  adjoint  Volterra  scattering 
operator  are  to  be  noted. 

We  have  seen  in  Chapter  I  that  for  the  absorption-only  case 

(a  -  y<)F(E,y;ic)  =  0  (2.7) 

the  spectrum  is  the  domain  C  in  which  (a  -  \ik)   vanishes.   The  singular 
eigenfunctions  were 


F(E,y;K)  =  F''"(K;E,y)  =  X(K)6(a  -  yK)  (2.8) 


where  the  second  identity  occurs  since  Eq.  (2.7),  the  streaming  equation, 
is  self-adjoint.  Thus  in  the  limit  of  no  scattering,  the  eigenfunctions 
of  the  fast  WBE  tend  to  the  delta- function  form,  Eq.  (2.8). 
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We  observe  from  Eqs.  (1.7)  and  (2.6)  that  the  domain  C,  due  to  the 
streaming  operator,  also  is  contained  in  the  continuous  spectrum  for  the 
slowing-down  WBE;  we  now  show  that  in  fact  it  is  identically  the  spectrum 
since  the  scattering  operator  will  cause  no  additional  contribution  to 
the  spectrum.  To  demonstrate  this  we  show  that  all  k  (6  C  are  in  the 
resolvent  set,  which  is  the  complement  of  the  spectriim,  and  is  defined 
as  those  values  of  k   for  which  (A  -  k)  has  a  bounded  inverse.  Therefore 
we  consider  the  existence  of  solutions  to  the  equation 

(A  -  K)4>   =  s(E,y).  (2.9) 

We  examine  first  the  case  of  isotropic  scattering,  for  which  the 
scattering  operator  becomes 

^1       /■<»  n       »00 

J  du'  J  dE'  i:g(E',y'  ^  E,ii)-  =  J  dM'  J  dE'  ^  E^(E^  ^  E)-  .      (2.10) 
-IE  -IE 

Using  this  operator  we  may  write  for  Eq.  (2.9)  the  equivalent  equation 


(a  -  yK)(j)(E,y) 


rl 

dy 


'  J  dE^  I  Z^(E'  -V  E)(})(E^y')  =  S(E,y).   (2.11) 


For  values  of  k  g  C  we  may  divide  by  (a  -  ]s<)   and  integrate  over  y: 
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I  (}>(E, 


y)dy  E  (|)(E) 


-1 


_dy_ 
a-yK  , 
E 


dE- 


2  ^s^^' 


E)(|)(E') 


dU 


S(E,y) 

a-yK 


(2.12) 


Since  k  j2  C,  both  integrals  over  y  exist  and  Eq.  (2.10)  is  of  the  form 


(()(E)  =  f(E) 


dE' 


2  ^s^^' 


^  E)(J)(E')  +  g(E)        (2.13) 


or 


1  -  f(E) 


dE'  ^  E  (E'-E) 
2   s 


())(E)  =  g(E).         (2.14) 


Provided  that  the  scattering  kernel  is  bounded,  the  Neumann  series 
inverse 


*(E)  =  I 


n=o 


.00 


f(E) 


dE'  -;r  E  (E'  ^  E)' 
^       s 


n 


g(E)        (2.15) 


always  exists  [58].  Thus  (A  -  k)  has  a  bounded  inverse,  and  we  have 
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the  result  that  the  complement  of  C  is  not  in  the  spectrum  of  A.  An 
identical  argument  applies  to  the  adjoint  operator. 

We  can  extend  this  result  to  anisotropic  scattering  by  making  a 
P  expansion  of  the  scattering  kernel  and  cj);  the  procedure  of  Eq.  (2.12) 
then  results  in  a  set  of  coupled  Volterra  equations  which  must  be 
inverted.  Thus  for  rather  general  scattering  kernels,  i.e.,  those  which 
can  be  developed  in  a  finite  bounded  P^^  expansion,  we  have  the  result 
that  the  spectrum  of  the  wave  Boltzmann  operator  consists  only  of  the 
continuum  C. 

Forward  and  Adjoint  Slowing-Down  Eigenfunctions 

Since  the  point  spectrum  for  the  slowing-down  problem  is  empty, 
there  will  be  no  regular  eigenfunctions  and  corresponding  space-  and 
E,y-separable  eigenmodes.  To  obtain  the  singular  eigenfunctions  corre- 
sponding to  the  continuous  spectrum  k  6  C,  we  may  apply  the  technique  of 
Chapter  I.  Adding  A(K)(a  -  yK)6(a  -  yic)  to  the  right-hand  side  of 
Eq.  (1.12)  and  dividing  by  (a  -  yK)  we  find 


1    f^    ^°° 
-1     E 


dE'  E  CE'.y'  ^  E,y)F(E,y;K)  + 


+  A(K)6(a  -  yK),     k  6  C 

(2.16) 
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or  equivalently. 


a- UK 


.00 


dy' 


dE'  Z^(E%y'  ^  E,u) 


F(E,y;K)  =  A(K)6(a  -  mk) . 

(2.17) 


At  this  point  it  is  necessary  to  proceed  more  formally.   It  has  been 
observed  in  Chapter  I  that  the  integral  of  the  factor  (a  -  yK)   over 
E,y  exists,  since  it  is  a  pole.  We  would  like  to  extend  the  Neumann 
series  inverse,  which  we  used  in  Eq.  (2.15),  to  Eq.  (2.17).  Accordingly 
we  write 


F(E,y;K)  =  A(k)  I 
n=o 


a-y< 


(•1    '°° 
dy^ 


dE'  Z  (E',y'  -^  E,y) 
s 


n 


6 (a  -  UK) 


■1 


=   X(K)  $   •   6(a  -  yK) 


Fsp(E,y;K) 


K  e  c 


(2.18) 


as  the  forward  slowing-down  eigenfunction.  The  formal  "Case's  Method" 
derivation  of  Eq.  (2.16)  must  be  verified  for  specific  scattering  kernel 
models  by  means  of  more  careful  arguments  such  as  those  used  in  substan- 
tiating Eq.  (1.32)  [31,33,34];  it  appears  that  this  will  succeed  for 
"well-behaved"  scattering  kernel  models.  For  k   in  the  continuous  spec- 
trum of  A  the  inverse  of  the  operator  (A  -  k)  exists  but  is  singular  [35], 
so  it  is  with  some  justification  that  we  write  the  second  form  of 
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Eq.  (2.18),  defining  the  formal  inverse  scattering  operator  $.     Further, 
the  Neumann  series  expansion  has  an  interesting  physical  interpretation 
in  terms  of  familiar  iterated  collision  integrals. 

To  see  this  we  first  recall  that  the  zero  scattering  cross-section 
eigenfunction  of  Eqs .  (1.21)  and  (2.8)  represent  neutrons  streaming  with 
eigen-energy  E  and  direction  y  ;  this  delta-function  distribution  is 
also  the  n  =  0  term  of  F(E,y;K).  The  second  term  is 


Ff'^(E,y;<)  .-l^E^(E^.y^^E,y)    E  <  E^ 


=  0  E  >  E         (2.19) 

K 


which  may  be  interpreted  as  the  distribution  resulting  from  one  down- 
scattering  interaction,  multiplied  by  the  transport  factor  (a  -  \ik) 
which  is  peaked  at  E  and  y  .  Similarly,  higher  terms  in  the  expansion 
may  be  interpreted  as  the  result  of  n  down-scattering  interactions,  so 
that  the  entire  eigenfunction  may  be  regarded  as  the  result  of  excita- 
tion by  neutron  waves  streaming  with  E  and  y  ,  along  with  an  associated 
down-scattered  contribution  excited  by  the  streaming  portion.  The 
eigenfunction  is  nonzero  only  for  E  and  below,  since  only  down- 
scattering  can  occur.   (This  deduction  from  Eq.  (2.18)  is  valid  whether 
the  Neumann  series  converges  for  E  <  E  or  not.)  We  see  that  the  eigen- 
function singularity  consists  of  a  delta-function  contribution  and  a 
pole  contribution  at  E  ;  a  similar  structure  occurs  in  the  thermal 
continuum  eigenfunction,  Eq.  (1.32).  Also  we  note  that  in  the  iterated 
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integrals  each  singularity  is  smoothed  by  integration,  and  that  the 
unintegrated  pole  can  be  factored  out  from  each  term  of  the  series,  so 
that  we  suspect  that  the  Neumann  series  inverse  will  indeed  converge  for 
rather  general  classes  of  scattering  kernels. 

The  adjoint  eigenfunctions  may  be  obtained  by  an  identical 
procedure;  Eq.  (2.6)  leads  to 


t  f  V 

F  (k;E,u)  =  A'(k)  I 


n=o 


a-yK 


du' 


^n 


dE'  E  (E,y 


E%y^)- 


6 (a  -  UK) 


=  ^"''(k)  3^*  6  (a  -  UK) 


K  e  C 


(2.20) 


or 


F^(k;E,u)  e  fJ^(k;E,u) 


(2.21) 


where  A  (k)  is  an  arbitrary  complex  constant.  The  form  of  the  adjoint 

t 
Volterra  operator  requires  that  F  is  identically  zero  for  E  <  E  ; 

again  the  delta- function  and  scattering-associated  term  with  pole 

singularity  at  E  =  E  occur.  The  properties  of  the  forward  and  adjoint 

eigenfunctions  may  be  summarized  by  the  rearranged  expressions 


F(E,u;k)  =  A(K) 


6 (a  -  y<) 


1        "  fl        f  K         Z    fE-,y'  *  E,p)      1" 


•1        E 


x23CE^,y,-E',y') 


E  <  E 

—      K 


"*"(!<:)    [fi 


?■*"(<;£, y)   =  A'Ck)    I6(a  -  yK)   + 


1        fE         E    fE.y  ^  E^y')      1 


1        r  (^        (^         S^^'^ 


n 


y< 


■1  E 


K 


X  E^CE'.y'  ->  E^,y^) 


E  >  E 

—      K 
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=  0 


E   >   E 


(2.22) 


E  <   E 


(2.23) 
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Discussion  of  the  Slowing-Down  Eigenfunctions 

An  interpretation  of  the  forward  eigenfunctions  in  terms  of 
iterated  collision  integrals  excited  by  monoenergetic  unidirectional 
(i.e.  y  =  u  )  streaming  neutrons  already  has  been  given.  We  proceed  by 
considering  their  biorthogonality  properties.  In  general,  due  to  the 
condition  expressed  by  Eq.  (2.5),  forward  and  adjoint  eigenfunction 
pairs  corresponding  to  different  eigenvalues  are  orthogonal  under  a 
y-weighed  scalar  product.   For  the  same  eigenvalue  <,  Eqs.  (2.22)  and 
(2.23)  clearly  show  that  the  scalar  product  will  not  vanish,  due  to  the 
coincident  delta-functions.   (This  product  of  delta- functions  of  two 
variables  requires  careful  interpretation  in  terms  of  the  theory  of 
generalized  analytic  functions  or  some  other  approach;  for  an  introduc- 
tion to  the  literature  on  this  aspect  of  the  singular  eigenfunction 
technique  see  McCormick  and  Kul'cer  [29].)  The  biorthogonality  proper- 
ties of  the  slowing-down  eigenfunctions  are  illustrated  schematically  in 
Figure  2.1  in  terms  of  the  energy  variable.  The  eigenfunctions  must  be 
orthogonal  for  overlapping  energy-distributions  as  well  as  in  the 
trivial  case  when  the  distributions  are  nonoverlapping  in  energy. 

It  is  interesting  to  consider  the  expansion  of  a  monoenergetic 
source  in  slowing-down  eigenfunctions.  This  is  schematically  repre- 
sented in  Figure  2.2.   We  see  from  the  first  two  sketches  that  such  a 
source  will  excite  not  only  continuum  modes  having  the  eigen- energy  E  , 
but  also  will  excite  to  some  extent  all  modes  with  lower  eigen-energies. 
As  is  apparent  in  the  third  sketch,  continuum  modes  with  higher  eigen- 
energies  will  not  be  excited. 
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Figure  2.1    Orthogonality  of  forward  and  adjoint 
slowing-down  eigenfunctions. 
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Figure  2.2    Excitation  of  slowing-down  eigenfunctions 
by  a  monoenergetic  source. 
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An  analysis  of  the  static  slowing-down  transport  equation  has  been 
performed  by  Maclnerney  [59]  for  constant  cross  sections  and  elastic 
scattering,  in  the  lethargy  variable.  By  performing  a  lethargy  Laplace 
transform  he  reduces  the  lethargy-dependent  problem  to  one-speed  trans- 
port form.  For  the  transformed  problem  (for  slowing  down  in  hydrogen) 
both  discrete  and  continuous  spectra  arise,  as  is  usual  in  the  one-speed 
problem  (see  standard  works  such  as  Case  and  Zweifel  [25]).  However 
due  to  inversion  of  the  lethargy  transform,  the  discrete  modal  contribu- 
tion fails  to  give  a  space-separable  solution  for  the  isotropic  space 
and  lethargy  Green's  function  (i.e.  a  source  6(u)6(x)).  This  is  in 
accord  with  our  result  that  a  monoenergetic  source  excites  a  continuous 
distribution  of  eigenfunctions.  Maclnerney  tentatively  attributes  his 
continuum  eigenfunctions  to  streaming  first-flight  source  neutrons; 
confirmation  of  this,  and  further  correlations  between  his  work  and  the 
present  "exact"  method  must  await  more  detailed  investigation. 

The  existence  of  a  discrete  mode  in  the  lethargy-transformed 
problem  raises  an  interesting  point  with  respect  to  implementation  of 
the  continuum  singular  eigenfunctions.  A  dispersion  function,  associ- 
ated with  both  discrete  and  continuum  modes,  was  seen  to  arise  naturally 
in  the  treatment  of  the  thermal  problem.  We  may  associate  the  disper- 
sion function  with  inversion  of  the  Fredholm  thermal ization  operator, 
since  in  the  slowing-down  case  only  the  Volterra  operator  is  present, 
and  no  such  dispersion  function  appears.  Physically  we  distinguish 
between  energy- regeneration  which  can  occur  through  up- scatter  in  the 
former  instance  and  energy  degradation  in  the  latter.   In  the  presence 
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of  energy- regenerative  mechanisms  we  find. the  potential  for 
establishment  of  E,y-space-separable  modes  (for  moderate  frequencies  and 
absorptions)  with  attenuation  length  longer  than  the  neutron  mean  free 
path.   For  the  slowing- down  problem,  with  such  mechanisms  absent  we  have 


?  =  Re  K  =     ^  (2.24) 

K 


SO  that  all  modes  are  attenuated  precisely  as  are  the  streaming-waves 
with  which  we  associate  them. 

However  it  is  well  known  that  the  neutrons  themselves  (e.g.  for 
neutron  pulses)  are  not  attenuated  in  this  manner,  even  though  no 
separable  mode  of  propagation  exists.  Evidently,  therefore,  we  are  not 
to  regard  a  continuum  mode  as  observable  or  capable  of  being  excited 
individually,  since  the  neutrons  which  would  constitute  such  a  wave 
certainly  would  not  be  attenuated  according  to  the  streaming  mean  free 
path.  This  is  further  evidenced  by  the  fact  that  a  monoenergetic 
unidirectional  source  excites  modes  having  lower  eigen-energies  as  well. 
Apparently  the  identification  of  an  individual  mode  with  streaming  and 
associated  scattered  neutrons  must  be  applied  with  some  caution,  although 
it  is  clear  that  actual  streaming  source  neutrons  are  represented  by  the 
delta- function  term  of  the  appropriate  eigenfunction.  We  must  conclude 
that  the  spatially  persistent  nonseparable  neutron  population  (as 
opposed  to  uncollided  neutrons)  excited  by  a  delta-function  source  is 
represented  by  constructively  interfering  continuum  eigenf unctions, 
where  this  constructive  interference  is  due  both  to  the  distributed 
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part  of  the  eigenfunction  excited  by  the  streaming,  and  to  eigenfunctions 
of  lower  eigen- energies.  Evidently  the  discrete  mode  in  Maclnemey's 
transformed  problem  corresponds  to  this  constructively  interfering  modal 
contribution. 

It  should  be  noted  that  the  idea  of  interference  of  neutron  waves 
is  not  new,  having  been  postulated  as  early  as  1964  on  the  basis  of 
diffusion  theory  by  Perez  et  al .  [60]  to  explain  phenomena  observed  in 
wave  experiments  in  subcritical  assemblies.  More  recently,  in  the 
transport  treatment  of  polycrystalline  materials  by  Yamagishi  [44], 
interference  effects  have  been  seen  to  arise  from  interaction  of  a 
continuiim  contribution,  due  to  neutrons  with  energies  below  the  Bragg 
cutoff,  with  the  higher  energy  neutron  population.   In  the  present  fast 
nonmultiplying  problem  we  have  seen  that  modal  interference  is  necessary 
to  describe  neutron  wave  propagation  in  all  but  purely  absorbing 
materials. 

In  the  scime  context  it  is  interesting  to  consider  elastic 
scattering  from  very  heavy  nuclei.   In  this  case  the  energy  loss  per 
collision  is  sufficiently  small  that  wave  propagation  in  such  a  medium 
is  essentially  monochromatic.  Thus  monoenergetic  analyses  may  be  per- 
formed such  as,  for  example,  those  of  Ohanian  et  al .  [61]  and  Paiano  and 
Paiano  [62] .   In  this  case,  due  to  the  energy-sustaining  model  of  the 
collision  process,  space-angle-separable  monoenergetic  eigenmodes  occur 
which  are  less  attenuated  than  Z  .  We  realize  that  in  the  actual  energy- 
dependent  problem  an  energy  loss  does  occur  with  each  scattering  interac- 
tion, so  that  only  continuum  modes  are  present;  nevertheless  these 


52 


continuum  modes  must  superimpose  in  such  a  way  as  to  yield  the  almost- 
separable  wave  behavior. 

The  macroscopically  elastic  scattering  kernel  model  of  the  above 
discussion  may  be  written 


E   (E^y^  ^  E,u)  =  E   (E')KCy'  ->  y)6(E  -  E') .        (2.25) 


This  kernel  also  is  noteworthy  because  it  is  not  bounded.   Clearly  our 
discussion  of  bounded  scattering  kernels  in  establishing  the  resolvent 
set,  K  jS  C,  does  not  apply  and  we  find  spectral  contributions  do  arise 
for  <  )S  C.   The  model  of  Eq.  (2.25)  is  discussed  in  Appendix  C,  along 
with  several  limiting  procedures  which  may  be  used  to  attempt  to  derive 
the  strictly  monoenergetic  case  as  the  limit  of  the  almost-monoenergetic 
case. 

We  conclude  the  discussion  here  by  observing  that  another  way  of 
viewing  the  problem  of  elastic  scattering  from  heavy  nuclei  is  to 
consider  a  detector  with  an  energy  window  AE  wide  enough  to  detect  all 
elastically  scattered  neutrons;  one  should  then  obtain  experimental 
results  which  are  in  accordance  with  monoenergetic  theory.  That  is, 
the  detector  response  should  show  an  asymptotic  exponential  signal  decay 
corresponding  to  the  momoenergetic  fundamental  mode;  this  detector 
response  is  the  physical  equivalent  of  solving  for  the  zeroth  moment  of 
the  flux  rather  than  the  flux  itself.   In  this  instance  we  must  agree 
with  Doming  and  Thurber  [40]  who  remark  in  another  context  that  in 
attempting  to  correlate  theory  and  experiment  one  can  be  mislead  by 
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considering  only  the  asymptotic  behavior  of  flux  solutions  rather  than 
their  moments . 

Realistic  Cross  Sections:  Nonmonotonic  vE 

In  Chapter  I  analysis  was  restricted  to  total  cross  sections  such 
that  vZ  is  monotonic.  This  was  done  because  the  continuum  values  of  < 
and  all  possible  E,vi  pairs  are  in  one-to-one  correspondence  for  monotonic 
vZ  ,  a  requirement  of  the  generalized  analytic  function  treatment  upon 
which  we  rely  for  completeness  results  in  the  separable  kernel  case. 
Here  we  explore  briefly  the  consequences  of  relaxing  the  monotonicity 
condition. 

For  this  case  degeneracy  of  the  continuum  results.  From  Eq.  (1.19) 
it  is  apparent  that  8(<)  will  assume  the  same  value  more  than  once  when 
vZ  is  not  monotonic.  This  is  illustrated  in  Figure  2.3,  where  it  is 
evident  that  for  the  same  value  of  6,  but  different  energies,  nondegene- 
rate,  singly  degenerate,  and  doubly  degenerate  regions  occur.  Higher 
degeneracies  may  result  from  more  rapidly  oscillating  cross  sections.  We 
exclude  the  case  of  constant  vE  ,  which  must  be  treated  separately. 

When  the  continuum  is  degenerate  the  coefficient  (a  -  yK)  in  the 
forward  and  adjoint  eigenvalue  equations  becomes  zero  for  more  than  one 
E,u  pair  at  each  degenerate  <  point.  Thus  in  Eq.  (2.16)  and  the  corre- 
sponding adjoint  expression  we  may  make  the  replacement 


M 

X(K)6(a  -  yK)  ->  y  A  (k)6  (a  -  yK) 
^   ^         '■^  m    m 
m=l 


M 

y  X   (K)6(E  -  E  )6(y  -  y  ).  (2.26) 

'^,     m  "^  ^     Km-^  ^     Km 
m=l 
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Figure  2.3    Degeneracy  of  the  continuum  due  to 
nonmonotonic  vl  . 
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Thus  from  Eqs.  (2.18)  and  (2.20)  we  have 


M 


F(E,y;K)  =  I     A^(K)  $-6^(a  -  yK)  (2.27) 

m=l 


and 


,  M   ,     , 

F"^(<;E,y)  =  I     A^(hc)  $^ -    6^(a  -  \xk)  (2.28) 

m=l 


for  an  M-degenerate  k.  Clearly  since  there  are  M  arbitrary  A's,  M 
linearly  independent  eigenfunctions  can  be  constructed.  An  obvious 
choice  is  to  set  the  A's  equal  to  zero  for  all  but  one  6  ;  we  define  the 
M  eigenfunctions 


^SD.m^^'^'^^  =  ^m^^^  ^  •  ^m^^  "  ^^^  (2-29) 


and 


fI^   ^(K;E,y)  =  xIm   J5^  •  6    (a  -   y<)  (2.30) 

bL),m         m         m 


which  we  notice  are  biorthogonal  when  the  forward  eigen-energy  is  less 
than  the  adjoint  eigen-energy  but  are  not  necessarily  biorthogonal 
otherwise.  Also  we  see  from  Eq.  (2.22)  that  the  forward  eigenfunctions 
will  have  pole  singularities  at  all  eigen-energies  E   less  than  the 

isTl 
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delta-function  eigen-energy  E  •  A  similar  structure  occurs  in  the 

^  "'  Km 

adj  oint  eigenf unctions . 


Fast  Multiplying  Media:   Zero  Scattering 
Cross  Section 

The  fast  multiplying  medium  problem  best  may  be  approached  by  first 
considering  Eq.  (1.13)  with  the  scattering  operator  absent.  Since  the 
fission  interaction  kernel  is  separable,  Eq.  (1.13)  then  becomes  identi- 
cal in  form  to  the  thermal  WBE  with  separable  kernel,  Eq.  (1.24),  which 
was  discussed  in  detail  in  the  first  chapter.   Identifying  x  with  M(E) 
and  vZ-(E')  with  Z  (E^),  we  may  write  down  immediately  the  results  for 
the  nonscattering  fast  multiplying  WBE  from  Eqs.  (1.29),  (1.30)  and 
(1.32).  Thus  we  find  that  the  discrete  eigenvalues  are  given  by  the 
dispersion  law 


A(K,oa)  -  I  - 


vE.,  — ^ 
f  a-yK 


=  0 


(2.31) 


and  the  corresponding  regular  eigenf unctions  are 


F(E,y;K.)  =  X(K.)  — ^^ 


J 


K  ^  C, 


(2.32) 


The  singular  continuum  eigenf unctions  are 


F(E,y;K)  =  A(K) 


X 


f   K 


a-yK 


X 


6 (a-yK) 


K  e  C 


(2.33) 
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Adjoint  eigenfunctions,  which  we  obtain  for  later  comparison,  are 
readily  found  to  be 


F  (K  ;E,y) 


,     vE  (E) 

A  (K.)  — ^ 

^  j'    a-yK. 


K  /g  C 


(2.34) 


and 


F'^(K;E,y)  =  A^(k) 


vZ^(E) 


a-UK 


ix(E,) 


A 


+  6 (a-yK) 


K  e  c 


(2.35) 


where  we  have  used  the  definition  of  x(E)  from  Eq.  (1.10).  We  note  that 
the  same  dispersion  function  occurs  in  both  forward  and  adjoint 
expressions . 

By  analogy  with  the  thermal  problem  we  expect  a  symmetric  pair  of 
eigenvalues  for  moderate  frequencies  and  absorption.  We  further  expect 
that  the  set  of  eigenfunctions  of  Eqs.  (2.32)  and  (2.33)  will  have  full 
and  half-range  completeness  properties  (although  strictly  speaking 
these  properties  were  demonstrated  for  a  symmetrized  kernel  in  the 
thermal  case;  a  similar  symmetrization  transformation  could  be  per- 
formed in  the  fast  case) . 


Discrete  Eigenfunctions  for  Fast 
Multiplying  Media 

We  now  turn  to  solution  of  the  fast  WBE  with  down- scattering,  as 

represented  by  Eq.  (1.13).   For  k  jg  C,  we  may  divide  by  (a  -  \xk)    and 


58 


invert  the  identity  minus  the  scattering  operator  (under  the  conditions 
which  were  discussed  previously)  to  obtain 


F(E,U;k)  =  3 


X 


a-yK 


-00 


dy^ 


dE'  vZ^(E')FCE',y';K)     (2.36) 


-1 


using  the  inverse  operator  defined  in  Eq.  (2.18).  Taking  the  scalar 
product  of  this  equation  with  vZ_,  we  find  that  the  condition  for 


solutions  to  exist  is  that 


A(k,(jo)  =  1 


vZ-,  $   •  -^^ 
f      a-yK 


(2.37) 


which  defines  the  dispersion  function  and  the  dispersion  law  for  the 
discrete  eigenvalues  k. .  The  expression  for  the  regular  eigenf unctions 
then  is 


F(E,y;K.)  =  A(K.)  3  •  —^^ 
3      J      a-yic. 


K  ^  C, 


(2.38) 


This  expression  may  be  compared  with  Eq.  (2.32);  making  use  of  the 
Neumann  series  interpretation  of  $,  we  see  that  the  presence  of  down- 
scattering  in  the  problem  has  resulted  in  an  addition  of  all  iterated 
collision  integrals  of  the  nonscattering  eigenfunction  (cf.  Eq.  (2.22)). 
Thus  the  discrete  eigenfunction  consists  of  the  fission  spectrum. 


59 


weighted  by  the  transport  factor  (a  -  uk)~^    (which  is  peaked  at  y  =  ±  l 
but  not  singular,  for  K  0  C) ,  and  smeared  down  in  energy  by  similarly 
weighted  scattering  operators.  We  will  discuss  the  regular  eigenfunc- 
tions  and  the  dispersion  law  in  more  detail  in  Chapter  IV. 

The  corresponding  adjoint  eigenfunctions  similarly  are  found  to  be 


t  ft    ^^ff^^ 

F  CK.;E,y)  =  A^Ck.)  $^   .  j^   .  (2.39) 


It  is  readily  verified  that  the  same  dispersion  law  is  obtained  here  as 
for  the  forward  problem. 

Continuum  Eigenfunctions  for  Fast 
Multiplying  Media 

By  application  to  Eq.  (1.13)  of  the  arguments  used  in  arriving  at 
Eq.  (2.18),  we  find 


F(E,y;K)  =  $ 


a-yK  j 
-1 


dy- 


dE^  vZ^(E')F(E',y^;K)   + 


+  A(k)   $  •  6 (a-yK)         <  e  C         (2.40) 
which  in  view  of  Eq.  (2.18)  may  be  written 

F(E,y;K)  =  $  .  -X_  f^j^^^f)  +  A(K)Fgp(E,y;K)   .      (2.41) 
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Eliminating  the  scalar  product  term  results  in  the  expression 


F(E,U;k)  =  A(K) 


X_  ^"^f'^SD^ 


a-VK 


+  F 


SD 


K  e  c 


(2.42) 


for  the  forward  continuum  singular  eigenfimction,  where  A  is  defined  by 
Eq.  (2.37).  The  adjoint  continuum  eigenfunction  is 


f'''(k;E,u)  =  a'^(k) 


vE, 


a-yK 


A 


+  F 


SD 


K  e  c 


(2.43) 


t   _   t 
where  in  this  case  F  =  $     •  6 (a  -  yK) . 


Discussion  of  the  Continuum  Eigenfunctions 

The  eigenfunctions  represented  by  Eq.  (2.42)  have  an  interesting 
interpretation  much  in  the  same  manner  as  that  of  Eq.  (2.22),  and  with 
similar  reservations  applicable.   Making  use  of  Eq.  (2.37)  for  the  dis- 
persion function  and  expanding  the  inverse  in  a  power  series,  Eq.  (2.42) 
may  be  written 


F(E,y;K) 

A(K) 


X 


a-yic 


1  + 


f     a-yK 


f    a-yK 


^^f'^SD 


-     F5p(E,y;<) 


(2.44) 
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Now  let  us  regard  (vE  ,F  )  as  the  initial  excitation  of  the  fission 
contribution  to  the  mode.  The  resulting  fission  neutrons,  after  being 
smeared  in  energy  by  scattering,  have  the  energy  and  angular  distribu- 
tion $   •  — ^  .  The  form  of  the  dispersion  function  expansion  suggests 

a-yK 

that  it  be  interpreted  as  a  modal  multiplication  due  to  the  sum  over  all 
generations  of  fission  neutrons .  Thus  we  see  that  the  continuum  mode 
again  may  be  regarded  as  streaming-associated,  since  it  consists  of  two 
terms  which  we  interpret  as  follows.  The  second  term  is  F   which  has 
been  seen  to  be  the  down- scattered  distribution  associated  with  stream- 
ing neutrons  having  the  eigen-energy  and  eigen-angle.  The  first  term  of 
Eq.  (2.44)  then  may  be  interpreted  as  the  fission-produced  modal  flux 
distribution  due  to  excitation  in  turn  by  the  scattered  term. 

This  attractive  exegesis  must  be  tempered,  as  in  the  slowing-down 
case,  by  considering  the  scalar  product  of  the  adjoint  eigenf unctions, 
Eq.  (2.43),  with  a  monoenergetic  source  function.  We  observe  first  that 
the  source  will  excite  modes  with  lower  eigen-energy,  due  to  the  term 

F  .   In  addition,  modes  having  eigen-energies  both  above  and  below  the 

^^  \)T. 

t  f 

source  energy  will  be  excited  due  to  the  fission  term  $  •  — -—  . 

"•^  a-yK 

These  results  may  be  compared  qualitatively  with  solutions  for  the 
static  fast  multiplying  medium  transport  problem  obtained  by  Nicolaenko 
and  Zweifel  [63]  and  Nicolaenko  [14].  Energy-transform  techniques  were 
used  to  treat  fission  and  elastic  scattering  with  constant  cross  sections 
in  the  former  study.   Inelastic  scattering,  the  model  for  which  already 
has  been  the  subject  of  comment  here,  was  added  in  the  latter.  Although 
detailed  comparison  again  is  difficult  due  to  the  complex  structure  of 
the  continuum  eigenf unctions,  we  find  consistencies  between  their  Green's 
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fimction  results  and  the  present  work.  Specifically,  in  both  studies. 
Green's  function  solutions  are  found  to  contain  both  space-separable 
contributions  (which  we  ascribe  to  the  discrete  eigenmode)  and  nonsepa- 
rable  "slowing-down  transients,"  which  are  solutions  to  the  slowing- 
down  equation  without  the  fission  term,  and  which  were  found  to  be 
necessary  to  achieve  completeness  for  the  eigenfunctions  of  the  appropri- 
ate Boltzmann  equation.  The  correlation  with  our  results  is  apparent. 

Degeneracy  of  the  Continuum 

Should  further  complications  seem  desirable  at  this  point, 
consideration  may  be  given  to  the  effect  of  degeneracies  in  the  continuum 
upon  the  above  treatment  of  continuum  eigenfunctions.  Since  the 
details  are  straightforward,  we  simply  note  that  linearly  independent 
sets  of  eigenfunctions  can  be  obtained;  in  particular  a  set  correspond- 
ing to  those  of  Eqs.  (2.29)  and  (2.30)  may  be  derived  by  an  identical 

procedure.  The  eigenfunctions  are  given  by  Eqs.  (2.42)  and  (2.43)  with 

t  t 

the  substitution  of  F_^   and  F     for  F   and  F  . 

SD,m      SD,m      bu  oU 

The  Boltzmann  Equation  with  Isotropic 
Interaction 

Finally,  some  general  consequences  of  isotropy  in  the  WEE  operators 
will  be  derived  for  use  in  Chapters  III  and  IV.  For  the  isotropic  ker- 
nel we  write 

K(E',u'  ^  E,y)   =  K(E'  ^  E)  (2.45) 


so  that  the  WBE  for  k  0  C  may  be  written 


F(E,y;K)    = 


a-yK  J 


f°° 


dy' 


dE'  K(E'  ^  E)F(E",y-';K:) 


■1  o 
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a-ytc  J 


KCE'  -V  E)F(E';ic)dE' 


(2.46) 


with  the  definition 


F(E;k) 


1 


F(E,y;K)dy 


(2.47) 


Integrating  Eq.    (2.46)   over  y  we  obtain 


F(E;k)      =     f(E) 


•*  N    1 1-.  "^ 


K(E     -y  E)F(E    ;K)dE 


(2.48) 


with  f (E)   defined  as 


f(E)     E 


dy 
a-yK 


-1 


(2.49) 
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Upon  solution  of  Eq.  (2.48)  we  then  reconstruct  the  angular  flux  from 


F(E,U;k)   =  f_I|lF(E;K) 


(2.50) 


For  K  e  C  the  continuum  eigenfunction  equation  becomes 


fl 


fCO 


F(E,u;k)  = 


a-y<  J 
-1 


du' 


dE'  K(E'  ^  E)F(E',y';K)  +  X(K)6(a-yK).    (2.51) 


Performing  integrations  over  y  we  have 


F(E;k)  =  f(E) 


K(E'  ^  E)F(E';K)dE'  +  X(k)6(E-E  ) 

1^ 


K  e  c, 


(2.52) 


CHAPTER  III 
APPLICATION  TO  THE  TRANSFER  MATRIX  METHOD 

Introduction 

In  this  chapter  the  analytic  results  obtained  for  the  WBE  will  be 
applied  to  the  transfer  matrix  formalism  of  Aronson  and  Yarmush  [64]  and 
Aronson  [65-70],  making  it  available  in  a  continuous- energy  transport 
representation.  There  are  two  aspects  of  this  technique  which  make  it 
attractive  as  a  potential  method  for  numerical  applications  of  transport 
theory.  First,  it  provides  a  convenient  general  framework  for  "problem 
solving"  in  terms  of  certain  basic  operators  (see  Aronson  [67]  for  a 
number  of  examples).  Second,  it  provides  an  explicit  method  for  obtain- 
ing transmission  and  reflection  operators.  As  we  will  see,  constructing 
some  of  the  required  operator  inverses  will  be  equivalent  to  determining 
the  half- range  orthogonality  properties  of  the  WBE  eigenfunctions. 
Although  it  is  not  possible  in  general  to  do  this  analytically,  numeri- 
cal inversion  techniques  certainly  may  be  employed,  so  that  the  transfer 
matrix  formalism  provides  a  straightforward  approach  to  this  difficult 
aspect  of  finite  medium  problems. 

The  transfer  matrix  for  slab  geometry  and  its  associated  eigen- 
value problem  are  derived  in  Appendix  A.  Essentially  what  one  must  do 
is  find  the  spectrum  and  eigenfunctions  of  a  certain  operator,  aS,   in 
whatever  representation  the  problem  is  formulated.  From  these 
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eigenf unctions  all  the  relevant  transfer  matrix  operators  may  be 
constructed,  as  well  as  transmission  and  reflection  operators.  Here  we 
will  obtain  the  a,  6,  and  08   operators  for  energy- dependent  wave  trans- 
port with  an  arbitrary  interaction  kernel,  and  then  show  that  for  iso- 
tropic scattering  the  eigenfunctions  of  a6   may  be  expressed  in  terms  of 
WBE  eigenfunctions.* 

Formal  Operator  Relationships 

The  operator  relationships  required  to  construct  the  transfer 
matrix  H  for  a  slab  of  width  x  may  be  summarized  as  follows: 


H  =  Se  ™  S  ^ 


(3.1) 


where 


"  B 

B 

+ 

'« 

1 

s  = 

2 

^ 

^ 

B 

B 

+ 

"c 

C 

+ 

0-1 

1 

s 

2 

C 

C 

+ 

(3.2) 


*Definitions  used  here  correspond  to  those  used  by  Aronson  in 
Ref.  [66]  and  earlier;  some  quantities  differ  by  a  factor  of  2  in 
Refs.  [69]  and  [70]. 
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-tA 


e      0 


+Tr 


(3.3) 


and  r  is  diagonal.  To  obtain  the  operators  B^,  C^  and  ?  one  must  first 
diagonalize  an  auxiliary  operator  o6: 


X'^   a6  X  =  r^ 


C3.4) 


n2 


where  T     is  diagonal;  then 


B^  =  X  ±  6  X  r"-^ 


(3.5) 


and 


c^  =  x"^  ±  r"^  x"^  a 


(3.6) 


The  explicit  wave  transport  representation  of  these  formal  relationships 
will  be  developed  in  the  following  sections. 

The  Operators  a6,  X,  and  X" 

The  operators  a  and  6  are  defined  as  the  sum  and  difference, 
respectively,  of  the  operators  a  and  3,  which  were  found  in  Appendix  A 
to  be 


«00 

fl 

a  - 

dE^ 

dy' 

, 

^ 

J  6(E'-E)6(y^-y) 


-Kg^(E^y'-E,y) 
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(3.7) 


0     o 


B  = 


dE' 


rl   r 

dy' 


-KggCE  ,-y  -E,y) 


(3.8) 


o     o 


where  K   is  the  kernel  of  the  Boltzmann  equation  interaction  operator 
(operator  A  of  Chapter  I),  including  scattering  and  fission.   In  an 
abbreviated  notation  Eqs.  (3.7)  and  (3.8)  become 


~   a   K 

a  = — 

y   y 


(3.9) 


^  =  ^ 


(3.10) 


so  that 


o  =  a  +  l  =  ---(K.*-K) 


(3.11) 


and 


6  =  a  -  3  =  -  -  -  (ic""  +  K") 

y  y 


(3.12) 


Then  their  product  is 
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05  =^-^  (k"   -   K") 


-  (a  -  K""  -  K") 


C3.13) 


in  terms  of  an  arbitrary  interaction  operator.  This  expression  is 
considerably  more  simple  when  K  =  K  ,  since  the  awkward  middle  term  of 
Eq.  (3.13)  vanishes.   In  particular  this  occurs  when  all  interaction 
processes  are  isotropic;  we  will  assume  this  to  be  the  case  throughout 
the  rest  of  this  chapter.  We  then  obtain 


;;7  _  a    2a  ~ 


dE' 


dp' 


^  6  (E'-E)  6  (y^-y)  -  ^K^^CE^  -  E) 

ly  y 


(3.14) 


which  is  the  form  we  will  consider  here. 

We  now  wish  to  obtain  the  operator  X  and  its  inverse  X   which  will 
diagonalize  o6   as  in  Eq.  (3.4).  This  may  be  done  by  first  finding  the 
spectrum  and  eigenfunctions  of  the  operator  a6 ,    and  then  constructing 
X;  X   will  be  constructed  in  a  similar  way  from  eigenfunctions  of  the 
adjoint  operator.  The  validity  of  the  diagonalization  will  of  course 
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require  that  the  sets  of  eigenfunctions  are  complete.  We  write  the 
eigenvalue  equation  for  06   as 


oS   '   X(E,u;Y^)  =  Y^X(E,y;Y^).  (3.15) 


Now  06  is  an  integral  operator  over  E  and  y,  and  we  will  find  that 
in  general  it  will  have  an  area  continuous  spectrum  as  well  as  a  possi- 
ble discrete  spectrum.  To  clarify  the  correspondence  between  the  eigen- 

2 
function  X(E,y;Y  )  and  the  operator  X,  let  us  consider  for  a  moment  the 

simpler  case  which  occurs  when  06   is  an  ordinary  N  x  N  matrix  (as  it  is, 

in  fact,  for  the  raultigroup  diffusion  representation) .  Then  its  spec- 

2  ~2 
trum  consists  of  the  N  discrete  eigenvalues  y.;  F  is  the  diagonal 

2 
array  of  the  y.,   and  X  is  the  corresponding  matrix  made  up  of  columns  of 

-2 
eigenvectors,  X..  -   X(E.;y.).  The  matrix  X  is  then  a  transformation 
ij      1  J 

from  the  basis  generated  by  the  eigenvectors  corresponding  to  the  indi- 

2  ~-l 

vidual  Y->  to  the  discrete-energy  space;  similarly,  X   is  the  inverse 

trans  format  i  on . 

In  the  present  transport  case,  the  situation  is  entirely  analogous, 

but  the  summations  over  the  discrete  spectriom  must  be  supplemented  by  an 

2 
integral  over  the  continuum  values  of  Y  ,   and  summation  over  the  E.  is 

~2 
replaced  by  integrals  over  E  and  y.  Thus  T     is  the  diagonal  operator 

consisting  of  both  the  discrete  eigenvalues  of  a&   (if  any)  and  the  con- 

2 
tinuous  spectrum,  y   .     The  operator  X  is  made  up  of  "columns"  of  eigen- 

2        2 
functions  X(E,y;Y  )  with  Y  as  the  "index";  it  will  involve  both  an 

integral  over  the  continuum  and  a  possible  sum  over  discrete  contributions. 
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In  other  words,  for  the  continuum,  X(E,ij;y  )  is  the  kernel  of  an  integral 

2 
operator  over  all  continuum  values  of  y   . 

Assuming  that  the  set  of  eigenfunctions  of  a&   is  complete,  X  may  be 

2 
regarded  as  a  transformation  from  the  basis  y  to  the  basis  E,u,  while 

the  operator  X"  is  the  inverse  transformation.  Writing  X  formally  as 


dy  X(E,y;Y  )  • 


(3.16) 


(the  integral  is  understood  to  include  the  sum  over  the  discrete 
spectrum,  if  any)  and  X   as 


x-1 


fl 


dE 


du  x"^(Y^;E,y) 


(3.17) 


the  left  and  right  inverse  relations  become 


X  ^X 


=  I  = 


dE 


dy  X"^(Y'^E,y)  J 


dY^  X(E,y;Y^) 


2    2   2 
dy  6(y  "Y^  ) 


(3.18) 


and 
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XX"  "^ 


I  = 


dY^  X(E',y';Y^) 


Y 


,.00 

dE 

0     0 


dy  x'-^(Y^E,u) 


fl 


dE 


du  6 (E-E') 6 (y-y") 


(3.19) 


2   ^2  2 

where  6 (y  -y'  )  is  either  a  Dirac  or  Kroneker  delta  for  y  in  the 

continuous  or  discrete  spectrum,  respectively.  The  first  of  these 

expressions  is  a  biorthogonality  relationship  for  the  two  functions 

2       -12 
X(E,y;Y  )  and  X  (Y  ;E,y),  while  the  second  is  a  closure  requirement  over 

y  e  (0,1]  (although  we  will  see  later  that  this  closure  relation  is 

essentially  a  full-range  condition) . 


Spectrum  and  Eigenfunctions  of  o6 

Using  the  explicit  expression  for  a6,   Eq.  (3.14),  we  may  write  the 
eigenvalue  equation,  Eq.  (3.15),  as 


Y 


^y 


?.] 

■                          ' 

a 

a 

' 

X  = 

77  "^  Y 

Y 

J 

[y       J 

[^ 

2a 
2 


dE' 


dy'  Kg^(E' 


E)X(E^y';Y  ) 


(3.20) 


Inspection  of  this  equation  shows  immediately  that  the  coefficients  of 
X  on  the  left-hand  side  will  give  rise  to  a  continuum  which  is  identical 
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to  the  domain  C  which  we  defined  in  connection  with  the  WBE;  that  is, 
when  Y  6  C  it  is  also  in  the  continuous  spectrum  of  06.     Further,  we  now 
may  apply  the  singular  eigenfunction  technique  to  obtain  an  expression 
for  the  continuum  X  eigenfunctions.  Noting  that 


f              ^ 

f 

s 

2 

a 

a 

\x 

—  +  Y 

■^ 

-  Y 

[y       J 

[M 

. 

2a 

I_.-l 


[a-uY   a+UYj 


-1 


(3.21) 


we  find  that  for  Y  6  C 


1      .       1 


a-yY      a+yY, 


1    f°°         f^ 
dE' 

o  o 


dy'  Kg^(E'  ^  E)X(E',y   ;y  ) 


A  (Y)  6 (a-yY) 


(3.22) 


where  we  explicitly  consider  only  Y*s  contained  in  the  half  of  C  which 
is  in  the  first  quadrant  of  the  spectral  plane;  including  the  other 
half  of  the  domain  C  gives  redudnant  results.   For  y  ^  C  we  have  simply 


X  = 


1 


a-yY   a+yYj 


.00 


dE'  I  dy'  Kgg(E' 


^  E)X(E',y^Y  )     (3.23) 


for  the  rest  of  the  spectrum  of  06.  Now  since  the  stipulation  has  been 
made  that  K_„  is  isotropic,  Eqs.  (3.22)  and  (3.23)  may  be  integrated 

DC 

immediately  over  y  to  give 
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X(E;y  )  =  £(E) 


dE'  K^pCE^  ^  E)X(E-;y^)  *  X(y)6(E-E^),     y  €  C 


BE' 


(3.24) 


X(E;y  )  =  f(E) 


dE'  Kgg(E'  ^  E)X(E';Y^),     y^C  (3.25) 


where  £(E)  is  the  same  function  defined  in  Eq.  (2.49).   But  these 
equations  are  the  eigenvalue  equations  for  the  WBE  with  an  isotropic 
kernel.  Thus  we  have  the  following  results: 

(a)  The  spectrum  of  the  isotropic  a6  operator  is  identical  to  the 
spectrum  of  the  related  Boltzmann  operator. 

(b)  The  y- integrated  eigenfunctions  of  the  two  operators  are 
identical. 

That  is,  identifying  Y  with  k. 


so 


c°°     fl 
dE' 

o    o 


X(E;y  )  =  F(E;k) 


dy'  Kgg(E'  ^  E)X(E',y';Y  ) 


(3.26) 


dE'  Kg^(E' 


E)X(E^Y  ) 


fCO 


dE^  K3g(E- 


E)F(E';k) 
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,00 


dE 


^  f'^^   Se 


(E'  -y  E)F(E^y^K) 


(3.27) 


Using  this  identity  in  Eqs.  (3.22)  and  (3.23)  and  comparing  with 

Eqs.  (2.46)  and  (2.52)  for  eigenf unctions  of  the  WBE,  we  find  that  for 

an  isotropic  K„_ 
BE 


X(E,y;<  )  =  F(E,y;K)  +  F(E,y;-<) 


(3.28) 


The  Inverse  Operator  X 


:;-! 


To  obtain  the  inverse  operator  X"  analytically  the  most  obvious 
approach  is  to  construct  it  from  the  eigenfunctions  of  the  operator 
adjoint  to  o6;    it  then  should  have  the  required  biorthogonality  proper- 
ties of  Eq.  (3.18).  Referring  to  Eq.  (3.19)  we  see  that  the  appropriate 
scalar  product  is 


f°° 

i^.^)  = 

dE 

dy(l>(E,y)i|;(E.y) 


(3.29) 


which  leads  to  the  adjoint  isotropic  eigenvalue  equation 


a     2 

'y 


rl 


X-\k   ;E,y)  =  2J  dE'  J  dy'  -5^  Kgg(E  -v  E')X"^  (K^E',y') 


o     o 


(3.30) 


Defining  the  function 
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(3,31) 


and  substituting  in  Eq.  (3.30), 


a     2 


v+r  2  „  .   2a 
X  (k  ;E,y)  =  -^ 

y 


dE' 


dy'  Kgg(E  -^  E')x'''(K^E',y') 


(3.32) 


where  ^^^(E  -»■  E')  is  the  kernel  of  the  adjoint  WBE  considered  in 
Chapter  II.  By  duplicating  the  development  of  Eq.  (3.20)  we  have 
immediately 


and 


X'^(K^E)  =  F^(k;E) 


x"''(K^;E,y)  =  F'''(K;E,y)  +  F"'"(-K;E,y) 


(3.33) 


(3.34) 


making  use  of  the  eigenf unctions  of  the  adjoint  WBE  from  Chapter  II. 
Then  X"-^(K^;E,y)  is  given  by  Eq.  (3.31). 


Diagonalization  Operators 

These  results  now  may  be  used  to  obtain  expressions  for  the 
operators  B^  and  C^  which  diagonalize  the  transfer  matrix  itself 
(Eqs.  (3.1)  -  (3.3)).   From  Eqs.  (3.5)  and  (3.6) 
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?;-! 


B^=X±6xr"EX±^ 


(3.35) 


and 


C^  =  X  ^  ±  f"-^  X  ^  a  E  x'-^  ±  ^  ^ 


(3.36) 


Since  F  is  diagonal,  elements  of  K   may  be  written 


C(E,y;K)  =  ^ 


fl 


dE' 


dy- 


[^Ci6(E--E)6(y^-y) 


o     0 


i  hE^^'  -  ^^ 


X(E',y;0 


y<  Y/-C   2- 

—  X(E,y;K  ) 


(3.37) 


making  use  of  Eqs .  (3.12)  and  (3.20).  Then 


B^(E,y;<)  =  ■^X(E,y;K^) 


(3.38) 


For  isotropic  K„„  we  note  that 
oc 


(a-y<)F(E,y;K)  =  (a+y<)F(E,y;-K) 


(3.39) 


so  that  from  Eq.  (3.28) 
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B+(E,y;K)  =  2F(E,y;±<). 


(3.40) 


Similarly  ^   is 


t^   =  [  dK^-4  6(<^-K) 


dE 
o   o 


dyx"-^(K"^;E,y) 


•00 

dE^ 

J   ; 

o    o 


1 


dy'-^^p-  6(E'-E)6(y^-y). 

y 

(3.41) 


so  that 


and 


r"^(K;E,y)  =  ;^x"^(<^E,y) 


C+(K;E,y) 


y  ,  ,   -„t   ^  2y  „t 
-^  (a±yK)X  =  ±  — ^  F 


aK 


(±<;E,y) 


(3.42) 


(3.43) 


making  use  of  Eqs .  (3.31)  and  (3.34),  and  Eq.  (3.39)  which  also  applies 
to  F  . 


Full-Range  Orthogonality  and  Completeness 

At  this  point  B^  and  C^  are  still  determined  only  to  within  a 
normalization  function  of  k.  The  normalization  relationship  may  be 
obtained  by  considering  the  operators  S  and  S   of  Eq.  (3.2).   Substi- 
tution of  Eqs.  (3.38)  and  (3.43)  for  B^  and  C^  into  the  left  inverse 
expressions  S  S  =  I  yields  the  appropriate  expression,  which,  with  a 
little  manipulation,  is  also  found  (as  for  the  monoenergetic  case  [68] ) 
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to  be  a  whole-range  biorthogonality  relation  for  the  WBE  eigenf unctions 
F  and  F  .  Similarly,  the  right  inverse  equation  (which  in  view  of 
Eq.  (3.1)  expresses  explicitly  the  physical  requirement  that  as  the  slab 
width  T  is  decreased  to  zero,  the  transfer  matrix  must  reduce  to  the 

identity)  SS"  =  I  is  found  to  be  a  whole-range  closure  expression  for 

+ 
F  and  F  .  Thus  we  see  that  the  validity  of  the  diagonalization  in 

Eqs.  (3.1)  -  (3.3)  depends  on  the  full-range  biorthogonality  and  com- 
pleteness of  the  WBE  eigenfunctions,  properties  which  may  be  established 
in  the  context  of  the  WBE  itself. 

Half-Range  Orthogonality  and  Completeness 

The  entire  transfer  matrix  H  now  has  been  obtained  without 
explicit  application  of  any  half-range  conditions.   (It  was  not  neces- 
sary to  solve  Eqs.  (3.18)  and  (3.19)  for  the  inverse  of  X.)  However, 
we  recall  that  the  transfer  matrix  relates  the  angular  flux  over 
ye[-l,l]  at  one  slab  surface  to  the  corresponding  flux  at  the  other.   It 
is  apparent  that  use  of  the  full  transfer  matrix  is  completely  equiva- 
lent to  applying  full-range  boundary  conditions  at  one  interface  to  an 
infinite-medium  eigenfunction  expansion  of  the  flux,  which,  of  course, 
completely  determines  the  flux  elsewhere  within  the  slab,  and  in  particu- 
lar at  the  opposite  surface.   In  situations  where  only  incident  fluxes 
are  known  at  interfaces,  the  transmission  and  reflection  operators  T  and 
R  are  more  appropriate  to  the  problem.   It  is  in  constructing  these 
operators  that  half-range  conditions  will  appear. 

The  transfer  matrix  formalism  provides  expressions  for  T  and  R  in 
terms  of  the  diagonalization  operators  B^  and  C^  (Eqs.  (A. 63)  through 
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(A. 66)),  all  of  which  are  defined  for  y  6  (0,1].  However,  inverses  of 

the  operators  B^  and  C^  occur  explicitly,  and  must  be  evaluated  for  the 

interval  (0,1]  if  one  is  to  compute  T  and  R.  Thus,  for  B  ,  we  have  the 
inverse  relationships,  using  Eq.  (3.40), 


'°°   f  1 
dE   dy  B;;^(K;E,y)F(E,u;±K')  =  6(k'-k)         (3.44) 

0    o 


and 


2 
2 


dK  F(E,y;±K)  B;^(<;E",y')  =  6 (E'-E) 6 (y"-y) .      (3.45) 


Solving  the  singular  integral  equation,  Eq.  (3.44),  for  B^  is 

equivalent  to  obtaining  the  weight  function  W(E,y)  and  normalization 

+ 
N(k)  for  half-range  biorthogonality  of  F  and  F  ,  since  if  they  are 

complete  on  the  half-range  we  identify 


2B;^(<;E,y)  =  N(K)W(E,y)F'''(±ic;E,y)  (3.46) 


so  that  Eq.  (3.45)  becomes  a  closure  relation  for  the  WBE  eigenf unctions. 

It  was  pointed  out  previously  that  for  most  realistic  interaction 
models  analytic  solutions  of  Eqs .  (3.44)  and  (3.45)  are  not  likely  to  be 
readily  forthcoming,  and  consequently  the  transmission  and  reflection 
operators  must  be  constructed  without  analytic  expressions  for  the 
required  inverses.   (In  this  respect  the  invariant  imbedding  method  is 
an  alternative,  as  it  provides  a  completely  different  formulation  for 
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R  and  T.  See  Pfeiffer  and  Shapiro  [71]  for  a  review  of  the  various 
approaches  to  transmission  and  reflection,  and  Mockel  [54]  for  an 
example  of  application  of  both  transfer  matrix  and  invarient  imbedding 
methods  to  thermal  neutron  wave  propagation.)  Nevertheless,  supposing 
that  some  approximate  numerical  representation  has  been  found  for  the 
WBE  eigenfunctions,  one  can  always  invert  the  diagonalization  operators 
numerically,  which  we  have  seen  to  be  exactly  equivalent  to  numerically 
determining  half-range  normalization  of  the  WBE  eigenfunctions.  After 
the  numerical  inversion  is  performed,  R  and  T  also  may  be  calculated 
numerically. 

Application  to  Fast  Neutron  Wave  Propagation 

We  have  seen  that  for  isotropic  interactions  all  the  essential 
operators  required  to  apply  the  transfer  matrix  technique  may  be  con- 
structed from  solutions  of  the  WBE  eigenfunction  problem.   Further, 
requirements  for  half-  and  full-range  closure  and  biorthogonality  of 
these  eigenfunctions  are  implicit  in  the  formalism.  These  results  and 
in  particular  Eqs.  (3.28),  (3.34),  (3.40),  and  (3.43)  extend  Aronson's 
static  monoenergetic  work  [65-70],  and  are  quite  general  with  respect 
to  the  energy-dependence  of  the  kernel. 

It  is  interesting  to  note  that  relationships  between  transfer 
matrix  and  Boltzmann  equation  eigenfunctions,  analogous  to  these  equa- 
tions, have  been  obtained  for  arbitrary  anisotropic  scattering  and  for 
azimuthally  dependent  problems  in  the  monoenergetic  case  [69].  While 
there  is  no  new  physics  per  se  in  the  transfer  matrix  approach  beyond 
that  contained  in  the  WBE,  and  thus  one  expects  a  direct  interrelationship 


82 


between  the  two  eigenvalue  problems  [68],  it  remains  to  be  seen  whether 
the  simple  form  of  Eqs.  (3.28)  and  (3.34)  will  hold  for  energy- dependent 
formulations  with  anisotropic  scattering.   It  should  be  stressed  that 
since  no  stipulations  other  than  scattering  isotropy  have  been  made  with 
respect  to  the  interaction  kernel  in  this  chapter,  the  results  obtained 
here  apply  equally  to  fast  and  thermal  regions,  with  or  without  fission, 
etc.  Furthermore,  the  form  of  the  term  a(E,(jo),  which  we  have  taken  to 
be  Z     +   — ,  does  not  enter  explicitly  and  therefore  generalization  is 
immediate  to  finite  transverse  dimensions  through  introduction  of  a 
transverse  buckling.  Of  course  the  static  case  w  =  0  is  included. 

Thus  the  equivalence  of  the  continuous  energy  transfer  matrix  and 
WBE  eigenfunction  approaches  for  isotropic  interactions  is  quite 
apparent,  so  that  we  may  regard  the  transfer  matrix  method  as  a  possibly 
convenient  framework  for  application  of  the  WBE  analysis.   In  this  sense 
we  now  have  in  the  results  of  this  chapter  a  complete  treatment  of  the 
analytic  basis  of  the  isotropic  scattering  energy-dependent  transfer 
matrix.  Application  to  fast  neutron  wave  and  pulse  propagation  thus  is 
a  matter  of  finding  suitable  means  of  implementing  the  formal  results 
obtained  in  Chapter  II,  using  the  basic  relationships  developed  here. 


CHAPTER  IV 


APPLICATION  TO  DISPERSION  LAW  AND  DISCRETE 
EIGENFUNCTION  CALCULATIONS 


Introduction 

Expressions  for  the  spectrum  and  eigenfunctions  of  the  fast 
multiplying  WBE  were  derived  in  Chapter  II.  Emphasis  there  was  placed 
on  application  of  the  singular  eigenfunction  technique  to  obtain  formal 
expressions  for  the  continuum  eigenfunctions.   In  this  chapter  we  will 
consider  the  mathematically  more  straightforward  topic  of  the  discrete 
spectrum  and  regular  eigenfunctions  of  the  wave  transport  operator. 

This  subject  now  is  quite  thoroughly  understood  in  principle  for 
the  separable  fission  kernel  model;  Travelli  [16]  has  presented  an 
essentially  complete  transport  multigroup  numerical  treatment  of  the 
fast  wave  slab-geometry  eigenvalue  problem  which  takes  into  account 
scattering  anisotropy,  through  a  P^   expansion,  and  delayed  neutrons. 
However,  apart  from  Travelli 's  work  and  the  present  investigation  [72], 
this  author  is  not  aware  of  any  other  numerical  fast  neutron  wave 
results  which  have  been  reported.   It  is  surprising  indeed,  in  view  of 
the  apparent  timeliness  of  fast  space-dependent  kinetics  studies,  that 
use  is  not  being  made  of  tools  such  as  these  to  investigate  neutron 
disturbance  propagation  in  detail. 

In  this  chapter  an  interrelationship  between  the  dispersion  function 
and  regular  eigenfunctions  is  made  explicit,  and  some  properties  of  the 
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dispersion  law  are  noted.  An  application  of  these  relationships  is 
made  to  the  dispersion  law  for  multiplying  media  with  isotropic  elastic 
slowing-down  and  inelastic  scattering.  Numerical  results  are  presented 
for  the  case  of  one  elastically  scattering  species. 

The  Dispersion  Function  and  Discrete  Eigenfunctions 

The  dispersion  law  concept  in  the  context  of  neutron  wave 
propagation  originally  was  introduced  by  Moore  [73]  as  a  relationship 
between  wave  frequency  and  complex  wave  length;  subsequently  this  con- 
cept was  rather  broadly  generalized  [8-10],  as  discussed  in  Chapter  I. 
We  have  seen  here,  in  an  exact  transport  treatment  of  the  fast  neutron 
wave  problem,  how  the  dispersion  function  occurs  in  the  structure  of 
solutions  in  multiplying  media,  as  indeed  it  does  for  all  energy- 
regenerative  media,  as  a  result  of  the  presence  of  a  Fredholm  integral 
operator.  We  now  will  proceed  to  formulate  the  discrete  mode  eigenvalue 
problem  for  a  separable  fission  kernel  in  a  way  which  both  appeals  to 
intuition  and  suggests  a  method  for  computing  the  dispersion  law. 

The  dispersion  functions  for  fast  multiplying  media  with  a 
separable  kernel  were  given  by  Eqs .  (2.31)  and  (2.37)  for  the  cases 
when  downscattering  is  absent  and  present  respectively.  These 
expressions  may  be  summarized  by 


A  =  1  -  (vZ^,  G)  (4.1) 


where 


G(E,p;k)  =  $  •  -^ —  (4.2) 
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with  $  defined  from  Eqs .  (2.17)  and  (2.18);  when  scattering  is  ignored 
3=1.  The  operator  $   in  Eq.  (4.2)  may  be  inverted  and  the  equation 
multiplied  by  (a-yK)  to  give 


1    f°° 
(a-UK)G(E,y;K)  -  |  dy'   dE"Z^(E^y'->E,y)G(E,y;K)  =  x  •   (4-3) 

-1     E 


Thus  the  function  G(E,y;ic)  is  the  solution  to  a  pure  slowing-down  WBE 

with  the  normalized  fission  spectrum  x  as  a  source,  corresponding 

itot  -<x 
(cf.  Eq.  (1.4))  to  the  space-  and  time-dependent  source  xe   e 

From  Eq.  (4.1)  we  see  that  K  is  an  eigenvalue  when  G(E,y;K)  satisfies 
(vZ-,  G)  =  1;  in  that  case,  comparing  Eq.  (4.2)  to  Eq.  (2.38),  we  see 
that  G  is  an  eigenfimction  of  the  total  WBE  with  the  particular  nor- 
malization A(k)  =  1.  The  term  A'-^  already  has  been  interpreted  in  the 
discussion  of  the  form  of  the  multiplying  medium  continuum  eigen- 
functions  as  a  modal  multiplication  factor;  here  we  may  regard  the 
dispersion  law  as  a  modal  prompt  criticality  condition  which  determines 
the  value  of  ic  for  which  a  solution  to  the  homogeneous  eigenfunction 
equation  will  exist. 

An  analogous  interpretation  of  the  dispersion  law  for  the  die-away 
experiment  in  fast  multiplying  systems  has  been  given  by  Moore  [8]  and 
Doming  [41],  generalizing  the  work  of  Storrer  and  Stievenart  [74]  who 
arrived  at  an  expression  for  the  fast  pulsed  neutron  dispersion  law  by 
considering  successive  generations  of  fission  neutrons.  The  particularly 
simple  forms  above  and  in  the  pulsed  neutron  case  for  the  dispersion 
function  and  eigenfunction  in  terms  of  vE^  and  x  are  a  direct  consequence 
of  the  separable  form  of  the  fission  kernel  [8,41,74].   We  will  see  how 
these  results  generalize  for  a  degenerate  fission  kernel. 
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Algorithms  for  Evaluating  the  Discrete  Spectrtim  and  Eigenfunctions 

Two  different  numerical  approaches  to  the  solution  of  the  eigenvalue 
problem  are  suggested  by  the  form  of  the  WBE  and  Eqs .  (4.1)- (4.3). 
Travelli  has  employed  both  [16,30],  having  derived  the  techniques  by 
means  of  computational  considerations.  The  first  method  is  a  direct 
numerical  solution  of  the  P^^  multigroup  representation  of  Eq.  (1.13)  for 
K  ^  C;  other  than  the  requirement  for  complex  arithmetic  this  approach 
is  straightforward  [30] .  An  alternative  procedure  is  suggested  by  Eq. 
(4.3).   In  form  it  is  a  familiar  slowing-down  equation,  for  which 
solution  techniques  are  well  established.  The  function  G(E,v;<)  is 
obtained  readily  for  a  particular  value  of  <   by  solution  of  Eq.  (4.3); 
the  dispersion  function  A  may  be  evaluated  by  Eq.  (4.1).   Zeroes  of  A 
for  a  particular  oo  then  may  be  found  by  application  of  a  complex  Newton- 
Raphson  procedure.  This  procedure  has  the  advantage  of  not  requiring 
solution  of  a  matrix  eigenvalue  problem,  which  can  become  prohibitively 
lengthy  for  the  large  numbers  of  energy  groups  needed  to  achieve 

accuracy  in  fast  medium  problems  [16] . 

3A 
The  Newton-Raphson  procedure  can  be  expedited  by  computing  r— 

a  K 

by  the  following  scheme  [73].  Differentiating  Eq.  (4.1)  yields 


|^=-(«vll>  t^-^) 


and  from  Eq.  (4.3) , 


(a-yK)  -^ 


du'l  dE'  E^(E',M'-^E,y)  |^  =  yG  (4.S) 


-1 
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which  is  a  slowing-down  equation  with  yG  as  a  source.  Solution  of  Eqs . 
(4.3)  and  (4.5)  can  be  carried  out  in  parallel. 

Extension  to  Degenerate  Kernels 

Only  a  slight  additional  effort  is  required  to  formulate  the 
eigenvalue  problem  for  a  degenerate  Fredholm  kernel  (e.g.,  multiple 
fissioning  species).  Following  our  development  for  the  separable  kernel 
model  we  obtain  instead  of  Eq.  (2.36) 


M       (m)    f     f        <•  X 

F(E.y;<)  =  J:  3  .  1—     dyM  dE'  vl/'"JF(E'.y';K) 
m=l     ^"^"^   i  i  * 


M 


^G(m)^y(m)^  F) 


m=l 


(4.6) 


using  obvious  notation  for  an  M-term  degenerate  fission  kernel, 
Reducing  this  to  the  matrix  equation 


M 


(y^"\  F)  =  I      (yC").  G^™))  (yf'"^  F) 
m=l 


(4.7) 


in  the  usual  way,  we  find  the  condition  for  existence  of  a  solution  is 
the  dispersion  law 


Det 


I  .  t(y(n),  G^™^)] 


=  A  =  0 


(4.8) 


where  the  quantity  in  brackets  is  the  M  x  M  matrix  having  elements 


(yCn)^  G^m)^  ^ 


(n)      y^"^ 
f   '     a-yK  J 


(4.9) 
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It  is  interesting  to  note  that  the  G  "^  are  solutions  to  M 
slowing-down  problems  like  Eq.  (4.3),  each  with  the  same  scattering 
operator  but  with  the  source  energy  distribution  x    characteristic 
of  the  m   species.   Eq.  (4.8)  is  more  complicated  than  previous 
expressions  for  the  dispersion  function,  but  evidently  we  may  retain 
our  interpretation  of  A   as  a  modal  multiplication  factor;  we  notice 
that  all  combinations  of  fission  of  the  n   species  due  to  down- 
scattered  neutrons  from  the  m   species  occur. 

Values  of  k  which  satisfy  Eq.  (4.8)  may  be  determined  by 
straightforward  extension  of  the  Newton- Raphson  scheme  discussed  above. 
The  eigenf unction  F(E,y;K)  then  may  be  reconstructed  by  means  of  Eq. 
(4.6);  the  coefficients  (y   ,  F)  are  the  elements  of  the  eigenvectors 
of  the  matrix  equation,  Eq.  (4.7),  and  the  functions  G^  "^  will  have  been 
evaluated  in  satisfying  the  dispersion  law.  Thus  we  have  achieved  a 
general  extension  of  the  separable  kernel  analysis  to  the  discrete 
spectrum  and  eigenfunctions  of  the  WBE  with  slowing-down  and  a  degenerate 
fission  (or  thermal ization)  kernel. 

As  a  postscript  to  the  above  discussion,  we  note  that  the  entire 
procedure  is  identical  in  the  case  of  the  adjoint  eigenfunctions,  with 
the  obvious  transpose  and  interchange  of  x    and  "^^n       ,    and  with  use 
of  $     rather  than  $.     Further,  it  is  easy  to  show  that 


(yf'^),  G^'"))  E  (vi:/"\  gW)  =  (G'^Cn)^  ^(m)^^        ^^^^^^ 
where 

t    t   "^^f 
G     =  $^   '   -^-  (4.11) 

■^    a-yK 
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so  that  should  one  wish  to  construct  forward  and  adjoint  solutions 
simultaneously  it  is  necessary  to  solve  Eq.  (4.8)  only  once. 

Isotropic  Elastic  and  Inelastic  Scattering 

The  slowing-down  equations  encountered  in  the  previous  sections  are 
solved  readily  by  any  of  a  number  of  methods  available  from  the  fast 
reactor  literature;  see,  for  example,  the  review  of  Okrent  et  al,  [75]. 
To  illustrate  the  method  described  above  we  will  develop  the  expressions 
for  a  continuous  slowing  down  [76,77]  model,  with  the  addition  of  a 
simple  inelastic  scattering  model  as  well. 

For  isotropic  scattering,  Eq.  (4.3)  may  be  written 


G(E;k)  -  f(E) 

E 


fOO 

1 


2lg(E'^E)G(E";K)dE'  =  f(E)x        (4.12) 


using  Eq.  (2.49)  for  f(E);  cf.  Eqs .  (2.48)  and  (2.52).   Differentiating 
with  respect  to  k,  we  obtain 


^!iri-f(B)fi.^CB-.B)faB-  =  =Iftfll       (^-3) 


9k 

E 


which  is  the  isotropic  equivalent  of  Eq.  (4.5).   (The  y-dependent  eigen- 
function  can  be  constructed  from  the  y-integrated  form  by  using 
Eq.  (2.50).) 

For  isotropic  inelastic  scattering  from  M  species  and  inelastic 
scattering  fron  N  levels  with  a  constant  energy  loss  AE  per  interaction 
we  have,  in  standard  notation,  the  interaction  operator 
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E/a 


lOO 

NE 

dE'  E    (E^^E)    • 

=      1 

m=l 

m 


dE'-ll 


Z   ^'"^(E^) 


E  (1-a  ) 


NI  '°° 

1 
n=l 


dE"  Z  .  '^"^(E')6(E-E'+AE  ) 
SI  n 


(4.14) 


Defining 


H(E;k) 


G(E;k) 
f(E) 


(4.15) 


and  using  Eq.  (4.14),  we  find 


E/a. 


H(E;k)  = 


NE 

I 
m  1 


m  1  y   (m)  ,^- 
2  ^se   ^^  ^ 


E-(l-a  )   f(E')H(E')dE'  . 


NI 

+  1     J  Z^.(E+AE^)f(E+AE  )H(E+AE  )  +  X 
n=l 


(4.16) 


and  a  similar  equation  for  (4.13).  Eq.  (4.16)  is  solved  most 

E 

conveniently  by  converting  to  the  lethargy  variable  u  =  In  -^,   and 

h 

differentiating  with  respect  to  u.   The  resulting  equation 


3H(u;k) 
3u 


NE 


9  ,1   ,  .  u.    1  "f 


m=l 


E   (™^(u) 
se   ^  -^ 

1-a 
m 


f(u)H(u) 


u 


u+lna 


m 


NI 

+  T  i  ^M-^  [Z  .^"^(u  )f(u  )H(u  )E(u  )]        (4.17) 
2  ^1  E(u  )  3u  "■  SI   ^  n"^  ^  n    n    n 
n=l   ^  n    n 
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is  readily  integrated  niimerically,  beginning  with  zero  lethargy  located 
above  the  maximum  fission  spectrum  energy,  e.g.,  10  MeV.  In  Eq.  (4.17) 
we  have  defined 


H(u;k)  =  "■',( —  .  (4.18) 


T^ 


Illustrative  Results:  Dispersion  Law  and  Eigenfunctions  for 
Single  Scattering  Species 

This  method  of  solution  was  applied  to  Eq.  (4.17)  and  a  related 

8G 
equation  for  -r— ,  for  a  one-term  elastic  scattering  kernel.  These 

expressions  were  implemented  as  explicit  difference  equations.  The 

Newton-Raphson  method  employing  Eqs .  (4.1)  and  (4.4)  was  used  to  find 

zeroes  of  the  dispersion  function.  All  cross  sections  were  taken  to 

be  constant.  Note  that  with  constant  total  cross  section  the  boundary 

of  the  continuum  does  not  change  with  frequency. 

Two  computed  dispersion  laws  are  shown  in  Figure  4.1,  illustrating 
the  effect  of  different  values  for  multiplication.  Values  of  k  and  co 
are  expressed  in  units  of  Z  .  Computations  were  based  on  scattering 
from  a  nuclide  of  mass  230.  A  Maxwellian  fission  distribution  having 
an  average  neutron  energy  of  1.98  MeV  was  used.  Here,  as  in  the  work 
of  Travelli  [16],  only  one  discrete  mode  per  frequency  was  found,  in 
accordance  with  the  no-scattering  case  (i.e.,  by  analogy  with  the 
thermal  separable  kernel  results).  As  would  be  expected  the  less 
strongly  multiplying  eigenfunctions  are  in  general  more  rapidly 
attenuated. 

Eigenfunctions  corresponding  to  individual  points  on  the  dispersion 
curves  may  be  reconstructed  from  the  computed  functions  H(u)  using 
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Figure  4.1    Dispersion  laws  for  constant  cross-section, 
elastic  scattering  model. 


93 


Eq.  (4.18).  With  arbitrary  normalization  constant  A  the  isotropic 
eigenfunction  F(u;<)  is 

F(u;k)  =  XG(u;k) 

=  AH(u;k)  f(u)  e'"  .  (4.19) 

Zero  frequency  eigenfunctions,  which  are  real-valued,  are  shown  in 

Figure  4.2,  corresponding  to  the  zero  frequency  points  in  Figure  4.1. 

The  eigenfunctions  are  normalized  to  coincide  with  the  fission  spectrum 

X(u)  at  lower  energies;  the  effect  of  down-scattering  is  noticeable  at 

higher  energies. 

Amplitudes  and  phases  for  the  complex-valued  eigenfunctions  are 

shown  in  Figures  4.3  and  4.4.  Amplitudes  have  been  plotted  with 

different  normalization  to  separate  the  curves.  These  and  subsequent 

data  refer  to  the  vZ.  =  0.09  dispersion  law.  For  relatively  moderate 

frequencies,  ■=—  <  1  x  10  ,  there  is  little  departure  from  the  zero- 

^t 
frequency  energy  distribution  for  energies  above  about  0.2  MeV.  With 

frequency  increasing  through  the  midrange  of  the  dispersion  law  it 


appears  that  neutrons  of  lower  energies  experience  increasing  difficulty 


keeping  pace  with  the  wave.  At  y~  =  4  x  10  ,  oscillations  in  the  phase 

't 


occur  at  about  0.5  MeV  accompanied  by  some  noticeable  structure  in  the 
amplitude.  As  frequency  increases  beyond  this  value  severe  phase 
lagging  occurs  for  lower  energy  neutrons;  Figure  4.5  shows  this  effect, 
with  phase  and  amplitude  values  superimposed.   It  is  interesting  that 
the  amplitude  maximum  occurs  in  a  region  of  rapidly  varying  phase.  As 
frequency  increases  toward  the  continuum,  the  eigenfunction  becomes 
increasingly  peaked,  as  shown  in  Figure  4.6.  This  bunching  of  the 
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neutrons  in  energy  occurs  at  an  energy  value  corresponding  to  the 
maximum  of  (a-K)   ,  so  that  we  begin  to  see  the  onset  of  continuum- like 
behavior.  The  phase  corresponding  to  Figure  4.6  is  not  plotted;  the 
eigenfunction  exhibits  lagging  of  phase  (which  will  result  in  rapid 
oscillations  of  phase  relationships  over  small  ranges  of  energy)  even 
more  severe  than  that  shown  in  Figure  4.5,  beginning  at  about  0.8  MeV 
and  extending  to  lower  energies. 

Peaking  of  the  eigenfunction  at  the  energy  which  maximizes  (a-ic)" 
can  be  understood  by  noting  the  y- dependent  form  of  the  eigenfunction, 
Eq.  (4.2).  Recalling  that  (a-yK)~  has  a  pole  for  y  =  1  at  the  edge  of 
the  continuum,  it  is  apparent  that  the  discrete  eigenfunction  will 
become  more  nearly  monoenergetic  and  forward-directed  as  cj  approaches 
the  critical  frequency,  due  to  the  presence  of  this  factor  (note  that  it 
is  also  a  coefficient  of  each  of  the  collision  integrals  in  $) . 

Evidently,  then,  as  o)  approaches  the  critical  frequency  it  will 
become  increasingly  difficult  to  excite  a  discrete  mode  without  strongly 
exciting  continuum  modes  as  well,  in  view  of  the  lack  of  phase  coherence 
of  the  eigenfunction  in  its  maximum  amplitude  region.  Recalling  pre- 
vious remarks  on  the  apparent  interference  of  continuum  eigenfunctions 
to  produce  almost-separable  mode  behavior,  it  is  probable  that  the  con- 
tinuum contribution  to  the  propagated  wave  will  be  important  in  this 
case  as  well.   It  would  be  interesting  if  a  high-frequency  wave  were  to 
assume  asymptotically  a  sharp  energy  distribution  as  in  Figure  4.6;  this 
does  not,  however,  seem  likely.  Further  investigation  of  this  question 
will  require  methods  for  determining  superimposed  contributions  of  both 
discrete  and  continuum  modes. 


CHAPTER  V 


SUMMARY  AND  CONCLUSIONS 


Summary 

This  dissertation  has  explored  the  transport  formulation  of  the 
fast  neutron  wave  and  pulse  propagation  problem  in  slab  geometry.  Both 
multiplying  and  nonmultiplying  media  were  considered.  The  spectrum  of 
the  wave  Boltzmann  operator  and  formal  expressions  for  discrete  and  con- 
tinuum eigenfunctions  were  obtained.   It  was  found  that  while  the 
presence  of  the  scattering  operator  prevents  direct  application  to  the 
fast  case  of  the  analytic  techniques  which  have  been  used  for  the  thermal 
regime,  a  consistent  extension  of  the  thermal  eigenfunction  results  is 
found  which  takes  scattering  into  consideration.   In  this  extension  an 
operator  appears  which  arises  from  inversion  of  the  identity  minus  the 
scattering  operator;  the  inverse  operator  has  a  Neumann  series  expansion 
for  eigenvalues  in  the  discrete  spectrum.  This  expansion  represents 
familiar  iterated  collision  intervals,  each  integral  weighted  by  a 
transport  factor  (a-yK)~  .  The  same  operator  inverse  appears  in  the 
continuum  eigenfunction  expressions;  while  the  situation  there  regard- 
ing the  existence  of  the  Neumann  series  inverse  is  more  tentative, 
apparently  there  is  physical  and  some  mathematical  justification  for 
interpreting  the  inverse  operator  in  this  way.  All  eigenfunction 
expressions  reduce  to  forms  consistent  with  more  well-established  thermal 
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results  in  the  limit  of  vanishing  scattering  cross  section.  Qualitative 
correspondence  was  noted  with  the  few  transport  studies  dealing  with  the 
static  transport  equation  in  fast  multiplying  and  nonmultiplying  media. 

A  basis  for  application  of  the  wave  transport  results  to  finite 
regions  was  established  by  generalizing  the  transfer  matrix  formalism  to 
energy-dependent  isotropic  interaction  models.  The  complete  equivalence 
of  the  isotropic  Boltzmann  equation  and  transport  transfer  matrix  eigen- 
function  problems  was  demonstrated,  and  specific  simple  relationships 
between  the  two  sets  of  eigenfunctions  were  obtained.  These  results  are 
more  general  than  the  scope  of  the  present  fast  wave  investigation; 
since  the  specific  form  of  the  interaction  kernel  does  not  enter  into 
the  expressions,  they  apply  equally  well  to  fast  or  thermal  regimes. 

Finally,  disperrion  law  and  discrete  eigenfunction  expressions  were 
considered  in  detail  for  separable  and  degenerate  fission  kernel  models, 
an  efficient  algorithm  for  fundamental  mode  computation  was  presented, 
and  specific  expressions  were  developed  for  an  isotropic  elastic 
scattering  and  inelastic  scattering  model.  Eigenfunctions  were  shown 
to  consist  of  a  superposition  of  solutions  to  wave  transport  slowing- 
down  equations,  each  having  the  fission  spectrum  for  one  fissionable 
species  as  a  source.  Eigenvalues  were  found  from  the  dispersion  law, 
which  was  seen  to  be  in  the  form  of  a  modal  criticality  condition. 
Illustrative  calculations  were  performed  for  a  one-term  isotropic  elastic 
scattering,  separable  fission  kernel  model.  The  discrete  eigenfunctions 
were  found  to  have  an  energy  distribution  very  similar  to  the  fission 
spectrum  at  low  frequencies,  with  small  phase  variation  over  the  total 
energy  range.  At  frequencies  approaching  the  critical  frequency,  the 
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eigenfunctions  exhibited  strong  peaking  in  energy  and  in  the  forward 
direction,  accompanied  by  rapid  phase  variation  with  energy  in  the 
vicinity  of  the  peak. 

Conclusions  and  Suggestions  for  Future  Work 

The  theory  of  fast  neutron  wave  and  pulse  propagation  has  received 
almost  no  attention,  and  the  application  of  transport  theory  to  the  fast 
neutron  regime  has  received  little.  While  the  difficulty  of  performing 
analysis  in  the  fast  regime  is  acknowledged,  nevertheless  this  situation 
is  both  surprising  and  unfortunate.   In  this  dissertation  an  attempt  has 
been  made  to  obtain  specific  results  where  possible,  and  to  establish 
directions  which  future  investigations  should  take. 

Theoretically,  the  situation  with  respect  to  the  discrete  spectrum 
and  eigenfunctions  appears  to  be  well  in  hand,  at  least  for  moderate 
frequencies.  This  is  not  to  say  that  interesting  work  does  not  remain; 
as  noted  before,  very  little  numerical  work  exists  on  wave  and  pulse 
propagation  in  fast  multiplying  media,  and  evidently  there  is  no  experi- 
mental data  available.  The  qualitative  effects  of  material  composition, 
cross- section  characteristics  and  other  parameters  upon  fundamental  mode 
propagation  have  yet  to  be  explored.   In  particular  it  should  be 
ascertained  whether  traditional  propagation  experiments  are  sensitive 
to  pareimeters  which  may  be  of  interest  in  fast  reactor  kinetics  and 
safety  analysis.  The  implications  of  the  peaked  energy  and  direction 
exhibited  by  the  discrete  eigenfunctions  for  very  high  frequencies 
should  be  developed.   It  is  possible  that  this  behavior  is  simply  an 
indication  that  discrete  modes  are  not  excited  strongly,  relative  to  the 
continuum  modes,  as  frequencies  approach  the  critical  frequency.  On 
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the  other  hand  if  waves  should  be  found  to  exhibit  such  behavior  for 
deep  penetration  this  would  be  an  interesting  result;  this  subject 
could  be  further  explored  theoretically,  and  studied  experimentally  by 
means  of  energy-sensitive  detectors. 

To  investigate  high  frequency  waves  in  multiplying  media, 
propagation  in  finite  regions,  and,  of  course,  propagation  in  nonmul- 
tiplying  media  it  will  be  necessary  to  take  into  account  the  continuum 
eigenfunctions  in  implementing  an  eigenfunction  expansion  approach  (or 
the  equivalent  transfer  matrix  formalism) .  Whether  such  an  approach  can 
succeed  in  practice  is  a  question  which  must  be  addressed  by  further 
investigations.  Either  of  two  generic  methods  which  suggest  themselves 
may  be  employed  in  treating  the  continuum;  one  may  attempt  to  extend 
analytic  techniques,  such  as  the  generalized  analytic  function  method, 
to  the  fast  case,  or  one  may  search  for  some  approximation  to  the  con- 
tinuum eigenfunctions  which  discretizes  the  continuum  and  which  is 
amenable  to  niimerical  implementation.  The  former  suggestion  involves 
significant  extrapolation  of  present  mathematical  techniques,  and  may 
well  prove  too  complicated  in  realistic  applications,  though  valuable 
for  phenomonological  understanding  through  treatment  of  modeled 
interactions.  The  latter  suggestion  has,  to  the  best  knowledge  of  this 
author,  no  precedent. 

These  comments  are  not  intended  as  a  criticism  of  the  singular 
eigenfunction  approach  as  applied  in  particular  to  the  fast  neutron 
regime,  since  the  method  is  far  from  fully  developed  for  thermal 
analysis.   Indeed,  it  is  not  yet  clear  whether  the  singular  eigenfunction 
method  will  prove  useful  as  a  framework  for  numerical  application  to 
practical  problems  [29,50]  in  the  thermal  regime  although  it  is  of 
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course  valuable  both  as  a  source  of  exact  solutions  against  which  to 
check  other  methods  and  as  a  means  of  understanding  transport  effects. 
It  is  apparent  from  the  forms  of  the  continuum  eigenfunctions  for 
both  fast  multiplying  and  nonmultiplying  cases  that  it  is  necessary  to 
understand  the  slowing-down  continuum  eigenfunctions  in  order  to  under- 
stand all  continuum  phenomena.  These  eigenfunctions  are  important  not 
only  as  the  appropriate  ones  for  the  important  problem  of  propagation 
in  passive  media,  but  also  since  they  evidently  are  associated  with  the 
modal  excitation  mechanism  for  the  continuum  for  multiplying  media.   In 
this  connection  it  should  be  remarked  that  the  interpretation  which  was 
given  for  continuum  modes  can  be  viewed  in  rather  general  terms  by 
recalling  that  the  "Case's  method"  eigenfunction  expression,  e.g.,  Eq. 
(2.16),  looks  very  much  like  an  inhoraogeneous  Boltzmann  equation  with 
a  delta-function  source j  interpretation  of  the  singular  eigenfunctions 
in  terms  of  modal  excitation  by  such  a  streaming  source  therefore  may 
be  illusory.  On  the  other  hand,  in  the  absorption-only  case  it  is 
obvious  that  this  interpretation  is  correct.  Whether  this  is  an 
accurate  or  helpful  approach  to  understanding  the  continuum  is  a  matter 
which  can  be  resolved  only  by  detailed  investigation,  just  as  the  formal 
eigenfunction  expressions  themselves  remain  to  be  substantiated,  either 
theoretically  or  by  application  and  comparison  with  other  methods. 
Since  slowing-down  eigenfunctions  appear  in  all  continuum  expressions, 
a  logical  starting  point  for  investigation  of  continuum  phenomena  would 
be  to  explore  the  slowing-down  problem  for  a  simple  elastic  scattering 
or  comparable  kernel.  Results  of  such  a  study  would  be  not  only  of 
theoretical  interest  but  also  helpful  in  interpreting  pulse  propagation 
experiments  in  heavy  scattering  materials  [11-13]. 
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Mention  should  be  made  of  the  possibility  of  extending  the  present 
analysis  to  systems  having  finite  transverse  dimensions,  and  to  other 
geometries.  The  essential  analytic  novelty  presented  here  is  the 
inversion  of  the  scattering  operator;  applying  this  technique  to  other 
transport  problems  (again,  see  McCormick  and  Kuscer  [29])  should  be 
straightforward.  Also  it  should  be  repeated  that  the  methods  presented 
here  are  applicable  to  thermalization  problems  in  which  separable  or 
degenerate  kernels  are  employed. 
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APPENDIX  A 

Introduction 

In  Chapter  III  various  transfer  matrix  operators  in  slab  geometry 
are  expressed  in  terms  of  WBE  eigenfunctions,  and  the  basic  equivalence 
between  the  eigenfunction  expansion  and  transfer  matrix  techniques  is 
emphasized.  Actually,  the  transfer  matrix  approach  is  a  family  of 
methods  which  may  be  applied  to  fairly  arbitrary  geometries  and  to  physi- 
cal quantities  other  than  the  neutron  flux.   In  this  appendix  the  trans- 
fer matrix  will  be  derived  from  basic  principles  and  its  general 
properties  will  be  discussed,  essentially  following  Aronson  and  Yarmush 
[64]  and  Aronson  [65] .   It  will  be  found  necessary  to  introduce  the 
geometry  explicitly  in  order  to  progress  beyond  fundamental  considera- 
tions. The  most  general  expression  for  the  basic  wave-transport  transfer 
matrix  operators  in  slab  geometry  then  will  be  obtained. 

General  Formalism 

The  transfer  matrix  approach  to  neutron  propagation  differs  from 
the  more  common  Boltzmann  equation  technique  in  that  the  latter  is  con- 
cerned with  a  pointwise  description  of  the  transport  process  within  a 
medium,  while  the  transfer  matrix  method  considers  the  problem  of  trans- 
port through  a  medium  in  terms  of  incident  and  emergent  fluxes  at  the 
boundaries.   In  this  respect  it  greatly  resembles  the  invariant  imbedding 
technique,  and  in  fact  there  is  considerable  common  ground  between  the 
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two  [71].   Initially  we  consider  a  region  o£  space  occupied  by  a 
homogeneous  isotropic  medium  through  which  a  neutron  flux  may  propagate. 
(Here  the  flux  will  be  the  usual  energy  and  direct ion- dependent  "angular 
flux"  in  either  time-dependent  or  Fourier- transformed  representation. 
Obviously,  more  restricted  models  also  may  be  treated.)  All  that  is 
required  at  this  point  is  that  one  be  able  to  assign  two  "sides"  to  the 
region.  These  sides  will  be  arbitrarily  referred  to  as  "left"  and 
"right."  No  problem  arises  with  the  formalism  in  the  treatment  of  voids 
and  reentrant  configurations,  provided  that  all  relevant  surfaces  of  the 
region  under  consideration  are  assigned  to  either  the  left  or  right 
sides.   In  some  cases,  as,  for  example,  for  finite  bodies,  designation 
of  the  two  "sides"  may  be  entirely  arbitrary,  but  for  such  bodies,  once 
boundaries  are  assigned  the  transfer  matrix  is  defined  unambiguously. 
Since  only  slab  geometry  will  be  considered  in  detail  here,  this  aspect 
of  the  method  will  not  be  pursued;  it  will  be  assumed  that  appropriate 
sides  may  be  designated  for  all  regions  of  interest  [78] . 

The  transfer  matrix  may  be  defined  as  the  linear  operator  which 
relates  entering  and  emerging  fluxes  at  one  side  of  a  region  to  the 
corresponding  fluxes  at  the  other  side.   This  situation  is  represented 
schematically  in  Figure  A.l.  Here,  for  example,  (f)  (L)  =  (j)  (r  ,E,$^,t), 
where  r  is  the  coordinate  of  a  point  on  the  left  surface  of  the  region; 
motion  in  the  positive  sense  is  taken  to  be  from  left  to  right.   Thus 
(fi  (L)  and  <f)_  (R)  are  entering  fluxes,  while  <p    (L)  and  (|)  (R)  are  emerging. 
Entering  and  emerging  may,  of  course,  be  defined  unambiguously  by 
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</>H.(L) 


<^_(L) 


<^+(R) 


<^JR) 


Figure  A.l    Entering  and  emerging  fluxes  for  a 
single  region. 


<^+(L) 


i>JL) 


f^(R) 


<^_(R) 


Figure  A. 2    Entering  and  emerging  fluxes  for 
adjacent  regions. 
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j^  9  j^  •  n  >  0    =^ 


Q  3  '^  '   n     <  0         => 

s 


emerging 
entering 


where  n  is  the  outward-directed  unit  vector  normal  to  the  surface  at 
any  point. 

The  transfer  matrix  equation*  for  Figure  A.l  may  be  written 


"  <}>^(R)  " 

=  "l-r 

■  <^^(l) ' 

(|)_CR) 

Lj      l\ 

*_(L) 

CA.l) 


or  more  compactly. 


(t>(R)  =  \^j^<j)(L) 


(A.  2) 


where  the  transfer  matrix  H  is  a  2  x  2  matrix  of  linear  operators  when 
the  fluxes  are  resolved  into  entering  and  emerging  fluxes,  as  in 
Eq.  (A-1).   Linearity  of  H  requires  that  no  flux- independent  sources  be 
imbedded  within  the  region,  although  fission  sources,  which  are  linear 
functions  of  the  flux,  may  be  taken  into  account  in  H.  This  linearity 
condition  has  further  implications  with  regard  to  the  nature  of  the 
region  described.  Since  nonfission  sources  must  be  exterior  to  the 


*In  this  appendix  symbols  for  matrix  and  operator  quantities  will 
be  written  without  a  tilde,  since  there  is  no  possibility  of  confusion 
between  an  operator  and  its  elements. 
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region,  it  is  necessary  that  the  region  be  well  defined  by  its  bounding 
surfaces.  These  topological  nuances  will  not  be  of  interest  here;  it  is 
sufficient  to  note  that  Eq.  (A.l)  certainly  will  be  applicable  to  a 
singly-connected  region  totally  bounded  by  its  "right"  and  "left"  sides. 

Algebra  of  the  H-Matrix 

Consider  two  adjacent  regions,  1  and  2,  completely  bounded  by 
"left"  and  "middle,"  and  "middle"  and  "right"  sides,  respectively. 
This  situation  is  diagrammed  in  Figure  A. 2.  Applying  Eq.  (A.l)  to 
regions  1  and  2  we  may  write 


<j)(M)  =  H^(})(L) 

(J)(R)  =  H24>(M)  .  (A.  3) 


Combining  these  two  equations. 


(j)(R)  =  H2H^(()(L)  .  (A. 4) 


Evidently  the  transfer  matrix  for  the  combined  regions  1  and  2,  H  ,  is 


given  by  the  matrix  operator  equation 


H^2  =  H2»l  ^^-'^ 
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and  similarly  * 

H  regions  =  H  H  ,  •  *  •  H„H,  .  (A. 6) 

n"       nn-1        21  ^       ^ 

Obviously  transfer  matrices  do  not  in  general  commute.  Furthermore,  it 
is  apparent  that  for  any  region,  in  the  limit  of  vanishing  width  the 
transfer  matrix  must  tend  to  the  identity  operator: 


lim    H  =  I.  (A. 7) 

"width"  -^   0 


These  characteristic  properties  must  be  satisfied  by  all  transfer 
matrices  [64] . 

Form  of  the  H-Matrix:   T  and  R  Operators 

We  will  now  express  the  transfer  matrix  in  terms  of  linear 
operators  which  relate  emerging  fluxes  to  incident  fluxes.  This  may  be 
done  without  any  further  restrictions  on  the  type  of  region  considered. 
Figures  A. 3  and  A, 4  illustrate  the  situations  relevant  to  the  ideas  of 
transmission  and  reflection,  respectively.  An  incident  flux  0  (L)  will 
give  rise  to  an  emergent  flux  ({)  (R)  which  may  be  regarded  as  transmitted, 
and  4)_  (L)  which  may  be  regarded  as  reflected.   In  other  words,  any 
entrant  flux  at  one  side  will  produce  an  emergent  flux  distribution  over 
the  entire  surface;  the  flux  emerging  from  the  same  side  is  termed 
"reflected,"  and  flux  emerging  through  the  other  side  is  "transmitted." 
These  distinctions  apply  even  though  the  surface  configuration  may  be 
such  that  the  terms  are  not  particularly  descriptive  of  the  physical 
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<^+(L) 


<^+(R) 


Figure  A. 3         Transmission. 


</>+(L) 


<^_(L) 


Figure  A. 4         Reflection, 
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processes.  Transmission  and  reflection  operators  may  be  defined  for 
(j)  (L)  incident,  as  in  Figures  A. 3  and  A. 4,  by  the  relations 


4  (R)  =  T4  (L)  (A. 8) 


(t)_(L)  =  R(()^(L)  (A. 9) 


where  T  and  R  may  be  used  to  treat  both  attenuating  and  multiplying 
media.   These  operators  are  in  fact  Green's  functions  for  the  energy 
and  direction-dependent  flux,  expressing  outward- directed  fluxes  at 
surface  points  of  the  region  in  terms  of  the  "source"  (\)   (L) ,  which  is 
completely  arbitrary.  That  is, 

T  =  7(^^,\,E^  -^  P^,n^^,E^^:t^-t)  (A.  10) 

and 

R  E  R(^^,ti^,E^   ^  r'^,$^'j^,E'j^;t'-t).  (2.11) 

A  corresponding  pair  of  operators  may  be  defined  for  an  incident  flux 
distribution  at  the  right,  (j)  (R) : 

(J)_(L)  =  T*(i)_CR)  (A.  12) 

(}>^(R)  =  R*(!)_(R)  (A.  13) 

where  the  asterisks  represent  transmission  and  reflection  in  the  sense 
of  right- to- left.  These  operators  in  general  will  differ  from  their 
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unstarred  counterparts.  Equations  (A. 8),  (A, 9),  (A, 12) ,  and  (A.  13)  may  be 
combined  to  give  the  emergent  fluxes  for  arbitrary  incident  flux  distri- 
butions over  both  faces : 


*  (R)  =  T({)  fL)  +  R*(t)  (R) 


(()_CL)  =  R4)^(L)  +  T*(|._(R) 


(A. 14) 
(A.  15) 


But  Eqs.  (A. 12)  and  (A. 13)  may  be  solved  formally  for  the  fluxes  at  the 
right  face  in  terms  of  those  at  the  left.  (In  the  following  it  will  be 
assumed  that  all  required  inverse  exist;  in  applications,  existence  may 
impose  restrictions  on  the  problem.  For  example,  subcriticality  may  be 
required.) 


({)_(R)  =  -T*"^R(t)^(L)  +  T*"^(J)_(L) 

<I>^(R)  =  (T-R*T*"-^R)(|)^(L)  +  R*T*"-^())  (L) 


(A. 16) 
(A. 17) 


preserving  the  correct  sequence  in  all  operator  manipulations. 

If  the  fluxes  are  now  written  in  vector  form  and  the  operators 
combined  into  a  matrix,  we  have 


<t)^(R) 

T-R*T*"-^R 

R*X*-1 

<J)^(L) 

(t)_(R) 

-T*"-^R 

T*-l 

(t._(L) 

(A. 18) 
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which  is  precisely  the  transfer  matrix  equation,  with  the  identification 


H  = 


T  -  R*T*  ^R     R*T*"-^ 


-T*  R 


(A. 19) 


This  is  the  most  general  form  of  the  transfer  matrix  in  terms  of  the 
four  transmission  and  reflection  operators  [64] . 

Two-Region  Transfer  Matrix 

Equation  (A. 5)  now  may  be  applied  to  the  transfer  matrices  for  two 
adjacent  regions,  H  and  H  ,  as  in  Figure  A. 2.  The  transfer  matrix  of 
the  combined  region,  H  „,  is  the  matrix  product  of  two  matrices  such  as 
in  Eq.  (A. 19);  the  elements  of  the  product  matrix  bear  the  same  rela- 
tionship to  the  combined  region  as  the  corresponding  elements  in 
Eq.  (A. 19)  bear  to  a  single  region.   For  example,  we  may  compute  the 
lower  right-hand  element 


12 


*_  1     *  *_ 1     *_ 1  *_i 


'-1 


*  *-l 

^2^ri 


(A. 20) 


Inverting  this  expression. 


12 


^2^ 


1  * 


(A. 21) 
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which  gives  the  operator  for  transmission  through  both  slabs  from  right 

to  left,  provided  that  the  required  inverse  exists.  Then,  if 

*  -1 
(I  -  R-jR-,)   may  be  expanded  in  a  Neumann  series,  we  have 


00 

*    *   _ 

T  =  T     y 

12   'l   ^ 
n=o 


( 


h\ 


"  T*  .  (A. 22) 


Physically,  Eq.  (A. 22)  may  be  interpreted  as  in  Figure  A, 5.  The 
total  transmitted  flux  is  the  sum  of  contributions  from  flux  which  has 
suffered  n  pairs  of  internal  reflections  before  exiting  through  the 
left  face.  This  is  a  well-known  result  which  has  been  obtained  previ- 
ously by  "particle-counting"  arguments  analogous  to  the  one  illustrated 
here. 

Internal  Sources 


Sources  external  to  the  region  of  interest  are  expressed  in  terms 
of  the  entering  fluxes  which  they  generate  at  the  surfaces  of  the 
region,  and  therefore  enter  implicitly  into  the  formalism  through  <^    (L) 
and  <^_  (R) .  Fission  sources  are  linear  functions  of  the  flux  and  there- 
fore may  be  included  directly  in  the  transmission  and  reflection  opera- 
tors. Thus  the  only  type  of  source  which  would  occur  as  an  inhomogeneous 
contribution  to  the  transfer  matrix  equation  would  be  a  flux-independent 
source  imbedded  in  the  region.  Such  internal  sources  may  be  taken  into 
account  by  expressing  the  sources  in  terms  of  the  emerging  fluxes  which 
they  generate  at  the  left  and  right  faces,  Q  and  Q  ,  respectively,  as 
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Figure  A. 5    Transmission  through  adjacent  regions, 


Q_^ 


^Q. 


Figure  A. 6    Internal  inhomogeneous  sources. 
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in  Figure  A. 6.   Equations  (A. 14)  and  (A. 15)  may  be  modified  to  include 
the  contributions  of  internal  sources  to  the  emerging  fluxes: 


<}>^CR)  =  T(t)^(L)  +  R*(()_(R)  +  Q^ 


(f._(R)  =  R(j)^(L)  +  T*4)_(R)  +  Q_ 


(A. 23) 
(A. 24) 


These  two  equations  may  be  combined  to  give  equations  in  cj)  (R)  and  (()  (R) 
corresponding  to  Eqs.  (A. 16)  and  (A. 17),  but  containing  the  source  fluxes 
Q  and  Q_ ;  the  resulting  equations  similarly  may  be  arranged  into  matrix 
form  as 


(1)(R)  =  H<t)(L)  +  Q 


(A. 25) 


where  the  inhomogeneous  source  vector  Q  is  found  to  be 


Q  = 


*  *  1 
Q.  -  R  T  ^  Q 


*-l 
R    Q 


(A. 26) 


Transfer  Matrix  for  Homogeneous  Slabs 

The  rest  of  this  development  will  apply  specifically  to  homogeneous 
slabs  of  finite  thickness  in  the  x  direction,  and  infinitely  extensive 
in  the  y  and  z  directions.   (The  condition  of  homogeneity  is  not  as 
restrictive  as  it  may  seem,  since  slabs  which  vary  in  composition  in  the 
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X  direction  may  be  approximated  arbitrarily  well  by  layers  of  homogeneous 
sub-slabs.)   It  is  assumed  that  all  incident  fluxes  and  all  source 
distributions  are 

(a)  constant  with  respect  to  y  and  z 

(b)  constant  with  respect  to  azimuthal  angle  about  x. 

In  this  case  the  angular  neutron  flux  becomes  a  function  only  of 
X,  E,  y,  and  t,  where  y  is  the  usual  cosine  of  angle  with  respect  to  the 
X  axis.  Values  of  y  will  be  taken  to  range  from  0  to  1  for  both  (})^  and 
(})  .  This  designation  is  convenient  as  well  as  unambiguous,  since  the 
information  regarding  the  "direction"  of  the  fluxes,  normally  contained 
in  the  sign  of  y,  is  transferred  to  the  subscripts  +  and  -. 
For  a  homogeneous  slab  with  these  symmetries  imposed 

T  =   T(y^,E^  ^  y-j^,E'j^;T;t'-t)  =  T*(yR,Ej^  -^  M\,E'^;r;t^-t) 

(A. 27) 

and 

R  E  R(yj^,E^  -y  y'^,E'^;T;t'-t)  =  R*(\.Er  "^  ^'r'^^r'^'^^"^^ 


(A. 28) 


where  x  is  a  slab  width.  Since  4)  and  (J)  do  not  depend  upon  position  on 
the  surface,  T,  R,  and  H  have  spatial  dependence  only  through  slab 
width,  as  expressed  in  the  above  equations.  Equality  of  the  forward  and 
backward  operators  is  the  result  both  of  the  symmetry  of  the  physical 
situation  and  of  the  stipulation  y  6  (0,1].  An  immediate  consequence  of 
this  is  that 
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H  =  H*  (A. 29) 


and  hence  for  two  adjacent  slabs,  1  and  2,  of  identical  composition 


»12  =  »12  =  ("2"l^*  =  »l"2  =  »1»2  t^-^O^ 


so  that  in  this  case 


"l2  =  V2  =  Vl  =  "21-  ^^-''^ 


This  equality  also  is  obvious  on  physical  grounds,  since  the  composite 
slabs  (1  +  2)  and  (2  +  1),  and  the  total  slab  (12),  all  are  indistin- 
guishable with  respect  to  neutron  transport  properties. 

The  commutation  property  of  Eq.  (A. 31)  is  useful  in  determining 
the  functional  dependence  of  H  upon  slab  width  for  homogeneous  slabs, 
when  width  is  regarded  as  a  variable  parameter.  For  slab  width  x  +  Ax, 

H(x  +  Ax)  =  H(x)H(Ax)  (A. 32) 

independent  of  the  manner  in  which  Ax,  which  we  now  may  regard  as  an 
arbitrary  increment,  is  selected.  The  variation  of  H  with  respect  to 
slab  width  then  may  be  written 

HCx  +  Ax)  -  H(x)  =  H(x)[H(Ax)  -  I]  (A. 33) 
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so 


HCx   +  Ax)    -  HCx) 

Ax  "  "'■''-' 


HCAx)    -   l" 

Ax 


(A. 34) 


We  now  follow  the  usual  procedure  of  passing  to  the  limit  Ax  ->■  0. 
Define 


W  =   lim 
Ax  ■>  0 


I  -  HCAx) 

Ax 


(A. 35) 


That  W  exists  is  at  least  plausible,  since  the  transfer  matrix  is 
required  by  Eq.  (A.?)  to  have  the  property 


lim  H(Ax)  =  I 
Ax  -^  0 


(A. 36) 


An  explicit  form  for  W  will  be  obtained  by  use  of  Eq.  (A. 35) 
In  the  limit  Ax  =  0,  Eqs.  (A. 34)  and  (A. 35)  give 


3H(x) 

3x 


WH 


(A.  37) 


(note  that  W  and  H  commute)  which,  with  the  initial  condition  H(0)  =  I 
yields  the  form 


H(x)  =  e 


■xW 


(A. 38) 
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for  a  homogeneous  slab  of  width  x.  W  is  an  operator,  whereas  slab 
width  now  appears  only  as  a  scalar  multiplier.  We  now  have  two  expres- 
sions for  H.   Equation  (A. 19)  becomes 


H  = 


T  -  RT"-^R     RT"'^ 


t'^^r 


-1 


(A.  39) 


taking  into  account  the  symmetry  of  the  homogeneous  slab  geometry 
(Eqs.  (A. 27)  and  (A. 28)).   Equations  (A. 38)  and  (A. 39)  may  be  exploited 
to  determine  an  explicit  expression  for  W.  Now  the  transmission  and 
reflection  operators  must  have  a  Taylor  series  development  in  terms  of 
slab  width  x;  in  the  limit  x  ->-  0  we  should  have  to  first  order  in  x 


xa  +  • • ' 


(A. 40) 


R  =  xB  + 


(A. 41) 


since  T  should  tend  to  the  identity  operator  and  R  should  tend  to  the 
null  operator  as  the  slab  width  becomes  vanishingly  small.  Once  again, 
a  and  3  are  operators  (which  we  must  deduce  from  other  considerations) , 
while  x  is  a  scalar.  When  these  expression  are  inserted  into  Eq.  (A. 39), 
the  corresponding  equation  for  H  becomes,  to  first  order  in  x. 


lim  H  = 
x  ->-  0 


xa 


-  x6 


x6 


I  +  xa 


(A. 42) 
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Then  from  Eq.  (A. 35) 


W 


a     -6 
3     -a 


(A.  43) 


It  should  be  noted  that  W  is  dependent  of  slab  width,  since  the 
operators  a  and  6  depend  only  on  the  properties  of  the  slab  medium. 
Thus  through  the  introduction  of  slab  geometry  it  is  possible  to  reduce 
the  transfer  matrix  from  a  form  containing  rather  formidable  Green's 
functions  in  all  space  and  velocity  variables  to  one  in  which  the  spatial 
coordinate  enters  only  through  slab  thickness,  and  only  as  a  scalar 
parameter. 

The  Operators  a  and  g 

The  transfer  matrix  equation  for  a  homogeneous  slab  extending  from 
0  to  X  may  be  written  compactly 


<^(x)  =  H(x)(})(0) 


(A. 44) 


Differentiating  and  making  use  of  Eq.  (A. 37), 


M^  =  M^HO)   =  -WH(x)*(0) 


(A. 45) 
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so  that  by  use  of  Eq.  (A. 44) 


^^  =  -  W4>(x)  .  (A. 46) 


This  equation  is  especially  interesting  since  it  is  a  differential 
operator  equation  for  the  flux  vector  at  the  right  surface  only.  While 
Eq.  (A. 46)  was  obtained  for  a  slab  of  width  x,  it  applies  equally  well 
to  the  flux  at  a  surface  located  at  coordinate  x,  embedded  within  a  slab 
of  arbitrary  width  x,  as  in  Figure  A. 7.  Thus  Eq.  (A. 46)  is  essentially 
a  pair  of  coupled  differential  operator  equations  for  the  forward  and 
backward-directed  angular  fluxes,  (|)  (x)  and  (^   (x) ,  at  x  within  the  slab. 
For  example,  the  equation  for  the  <p     component  is 


34'Jx) 


+ 


a*  (x)  +  6(t)  (x)  (A. 47) 


3x      "^+ 

It  is  obvious  that  the  operators  a  and  3  must  contain  the  physics  of 
neutron  interactions  in  the  slab,  and  hence  contain  cross  section, 
scattering,  and  multiplication  operators.  The  fact  that  Eq.  (A. 47)  is  a 
partial  differential  equation  at  a  single  interior  coordinate  x  enables 
us  to  determine  the  form  of  a  and  3  from  the  equivalent  differential  or 
differentio-integral  representation  of  the  problem. 

Once  a  and  g  have  been  determined,  Eqs.  (A. 38)  and  (A. 43)  give  the 
operator  matrices  W  and  H,  and  in  principle  these  could  be  applied  to 
computations  without  further  development.  However,  the  operator  W 
appears  exponentially  in  H,  so  that  the  utility  of  this  direct  approach 
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Figure  A. 7    Fluxes  at  an  internal  coordinate  surface. 
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depends  on  rapid  convergence  of  terms  of  the  series  (tW)  .  While  this 
may  be  practical  for  some  problems  in  which  the  full  transfer  matrix  is 
required,  especially  for  thin  slabs,  an  alternative  approach  is  to 
diagonalize  W. 

Diagonalization  of  the  Transfer  Matrix 

Suppose  a  basis  has  been  found  upon  which  W  is  diagonal;  that  is, 
for  some  matrix  S, 


S  ^WS  =  A  (A. 48) 


where  A  is  the  diagonal  operator  consisting  of  the  spectrum  of  W.  Then 


W  =  SAS"-^  (A.  49) 

and 

„     -TW     -tSAS"'^  ..  cr>^ 

H  =  e    =  e  (A. 50) 

so  that 

-tA  -1 
H  =  Se   S   .  (A. 51) 


That  is,  if  the  spectrum  and  the  eigenfunctions  of  W  can  be  determined 

(the  eigenfunctions  constitute  the  matrix  S) ,  the  exponential  behavior 

-tA 
of  H  with  respect  to  x  will  enter  only  through  e   ,  which  is  itself 

diagonal. 

Rather  than  attempting  to  solve  an  eigenvalue  equation  for  a  2  x  2 

matrix  of  operators,  it  is  more  convenient  to  consider  first  an  auxiliary 

eigenvalue  problem.  Define 
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a  =  a  +  6 
6  E  a  -  e 


(A. 52) 
(A. 53) 


and  their  product  operator,  a6.  Suppose  that  this  operator  may  be 
diagonalized. 


r^  =  x"-^  a6X  , 


(A.  54) 


where  T     is  the  diagonal  operator  consisting  of  eigenvalues  of  o6 .     Then 
it  may  be  verified  by  substitution  that,  formally  at  least,  the  follow- 
ing equations  result: 


A  = 


-r 


(A. 55) 


-i 


(A. 56) 


,-1 


(A. 57) 


where 
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B^  =  X   ±  C 


(A. 58) 


and 


Explicitly, 


C^  =  x"^  ±  C^ 


5  E  6xr"-^  =  a'hv 


(A. 59) 


(A. 60) 


W  =  SAS 


-1 


1_ 
4 


and 


H  =  Se-"^S-^  =  i 
4 


0   -r 


e-^^   0 


0    e 


tF 


c    c 

+ 


(A. 61) 


C      C 
+ 


.  (A. 62) 


These  are  the  diagonalization  operator  equations  which  are  discussed  in 
detail  in  Chapter  III. 


Transmission  and  Reflection  Operators 

In  most  applications  it  will  be  necessary  to  compute  transmission 
and  reflection  operators  as  well  as  (or  in  place  of)  the  full  transfer 
matrix  H.   In  fact,  the  most  utilitarian  feature  of  the  formalism,  aside 
from  its  generality,  may  be  that  it  provides  relationships  between  trans- 
mission and  reflection  operators  and  eigenfunctions  of  the  more  familiar 
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neutron  transport  equations.  Expressions  for  these  operators  in  terms 
of  the  diagonalization  operators  may  be  deduced  by  comparing  the  form  of 
H  as  a  function  of  T  and  R,  Eq.  (A. 39),  with  Eq.  (A. 62).  With  a  little 
manipulation  one  obtains 


T(T)  =  4C:^  e""^^  E'^Ct)  (A. 63) 


+       + 


and 


defining 


and 


R(T)  =  E  (t)  E^^(t)  CA.64) 


E+(T)  =  B^  +  B-  e""^^  U  e"^^  (A. 65) 


U  =  C  C^^  =  -  B^^  B  (A. 66) 


which  permit  computation  of  T  and  R  from  operators  related  to  the 
physical  properties  of  the  slab  medium  through  o6   and  the  basic 
operators  a  and  B  . 

Wave  Transport  Form  of  a  and  B 

To  obtain  a  and  3  in  a  wave  transport  representation,  we  return  to 
Eq,  (A. 47).  Writing  a  and  3  as  integral  operators  (recalling  that  T  and 
R  and  thus,  from  Eq.  (A. 40),  a  and  3  are  defined  over  y  6  (0,1])  with  ker- 
nels K  and  K^,  Eq.  (A. 47)  becomes 
oc      p 
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a(()^(x,y,E,a)) 

"9r 


dE' 


dy'  K^(E  ,y' 


E,y;(jo)(j)^(x,y  ,E  ,03) 


r°°    fl 


dE'  I  dy"  Kg(E  ,y"  ^  E,y;oa)(t)_(x,y",E  ,0)) 


(A. 67) 


Now  in  similar  notation  the  WBE  is 


y  -^  +  a(E,a3) 


(})(x,y,E,a))  =   dE' 
o    -1 


dy'  Kgg(E',y''  ->  E,y)(t)(x,y',E'',a)) 


(A. 68) 


including  all  interaction  terms  in  the  kernel  K„„.  Rearranging, 

DC 


If  =  -  J  (|)  +  ^  J  dE'  j  dy'  Kgg(E',y'  ^  E,y)(l)(x,y',E',aj) 


(A. 69) 


When  y  is  restricted  to  the  interval  (0,1],  <t)i]i)    =  ({>  (y)  .   Furthermore, 
we  have  defined  <f>  (y)  E  (j)(-y).   If  the  integral  in  Eq.  (A. 69)  is  divided 
into  separate  integrals  over  positive  and  negative  y',  and  if  the 
substitution  y"  =  -y'  is  made  in  the  negative  integral. 


3(j)^(x,y,E,(i)) 


—  ({)     +  — 

tOO 

dE' 

y    +      y 

dy'  Kgg(E',y' 


E,y)(j)^(x,y  ,E  ,00) 


00    n 


+  -   dE' 

y 


dy"  Kgg(E',-y"  ->  E,y)<|)_(x,y",E',a3) 


o     o 


(A. 70) 


132 


where  y,  y^,  y"  6  (0,1].  Comparing  Eqs.  (A. 67)  and  (A. 70),  we  find 
that 


a  = 


J 
o 


dE-  J  dy'  pL^  6(E'-E)6(y"-y)  -  \  Kgg(E^y'  -  E, 


y) 


and 


(A.  71) 


3  S 


dE' 


dy' 


yl^BE^^ 


',-y'  ^  E,y)l 


(A. 72) 


A  similar  treatment  of  the  equation  for  ^   (x)  gives  identical  results. 

The  kernel  K   may,  of  course,  contain  any  useful  interaction  model. 
Be 

This  extends  Aronson's  work  [64-70]  to  the  most  general  wave  transport 
representation  of  the  one-dimensional  slab  transfer  matrix. 


APPENDIX  B 


Singularity  of  Inelastic  Scattering  Kernel  Models 

In  this  appendix  an  inelastic  scattering  kernel  model  [79,80]  will 
be  discussed  which  gives  rise  to  an  additional  spectral  contribution  in 
the  adjoint  problem.  This  observation  is  due  to  Nicolaenko  [14],  who 
made  use  of  the  kernel  in  treating  fast  static  transport  problems. 

We  recall  that  in  general  the  total  isotropic  scattering  cross 
section  is  defined  by 


13(E)  = 


E  (E->-E')dE'  .  (B.l) 

s^ 


For  downscattering  (either  elastic  or  inelastic)  this  becomes 


E^(E)  =   Z^(E-E')dE'  .  CB.2) 


We  may  write  the  kernel  in  question  as 

K.^CE-^EO  =  g(E)  h(E')    E  >  E' 

=0  E  <  E"  CB.3) 

For  constant  total  inelastic  scattering  cross  section,  which 
Nikolaenko  stipulates,  we  may  write 
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Z  . (E)  =  const. 

SI 


=  E 


SI 


g(E)  H(E')dE- 


(B.4) 


so  that  the  constant  cross-section  condition  requires  that 


g(E)  h(E')dE'  =  1 


or 


J" 


(E  )dE  = 


g(E) 


(B.5) 


If  h(E)  is  bounded  as  E^-0,  g(E)  -^0  in  this  limit.  Thus  while  the 
forward  inelastic  scattering  operator 


dE'  E  .  (E'-^E) 

SI 


Z^.  h(E)  I  dE'  gCE^) 
E 


(B.6) 


has  a  bounded  kernel  for  E  >  0,  the  adjoint  inelastic  scattering  operator 


.E 


dE"  Z  . (E^E')  •  =  Z  .  g(E) 

51  SI  "^ 


dE  h(E') 


(B.7) 


clearly  does  not,  as  g(E)  is  singular.  Thus  our  inversion  of  the 
identity  plus  scattering  operator  fails  for  the  adjoint  WBE  employing 
this  model.   In  fact  operator  (B.7)  has  eigenfunctions  [14]  g(E) 
with  corresponding  eigenvalues  y  foi"  ^H  ^  with  ReX2.1;  hence  the  spectra 
of  forward  and  adjoint  Boltzmann  operators  evidently  are  different,  and 
the  entire  approach  based  on  construction  of  biorthogonal  eigenfunction 
sets  is  unsuccessful  in  this  case. 
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The  question  therefore  becomes,  is  this  unfortunate  property  a 
characteristic  of  inelastic  scattering  in  general,  or  is  it  an  artifact 
introduced  by  the  specific  model,  Eq.  (B.3)?  Fortunately,  the  latter 
is  true.  From  Eqs.  (B.4)  and  (B.5)  it  is  evident  that  the  stipulation 
of  constant  total  cross  section  at  E^,  together  with  the  form  of  the 
kernel  (B.3),  forces  g(E)  to  be  singular.  That  is,  it  forces  a 
divergent  differential  cross  section,  a  result  which  is  both  incon- 
venient, as  we  have  seen,  and  physically  incorrect,  since  inelastic 
scattering  does  not  occur  at  all  for  sufficiently  low  energies. 
(Nicolaenko  avoids  the  difficulty  of  nonphysical  results  due  to  this 
divergence,  as  he  does  not  apply  the  kernel  below  the  inelastic  scatter- 
ing threshold;  nevertheless  the  singularity  does  emerge,  as  we  have  seen, 
in  the  adjoint  problem.)  As  we  are  dealing  with  fast  systems  we  are  in 
fact  indifferent  as  to  the  precise  form  the  inelastic  scattering  cross- 
section  model  at  thermal  energies  and  below,  but  we  certainly  are  at 
liberty  to  allow  the  inelastic  scattering  cross  section  2^^(E)  to  vanish 
as  E^,  for  example,  in  some  neighborhood  of  E  =  0.  This  additional 
degree  of  freedom  entirely  resolves  the  present  difficulty;  retaining 
the  kernel  (B.3),  we  see  that 


2„. (E) 
gCE) 


SI  (B.8) 


0 

so  that  for  E  ^  0  all  that  is  required  is  that  E^^(E)  vanish  sufficiently 

rapidly. 

It  is  apparent  that  the  singularity  of  the  kernel  adjoint  to  (B.3) 
is  not  a  property  which  generalizes  to  more  accurate  models  of  inelastic 
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scattering,  and  therefore  the  behavior  of  scattering  kernels  at  the 
zero-energy  limit  should  not  necessarily  hinder  the  development  of 
transport  theory  for  fast  systems. 


APPENDIX  C 


MACROSCOPICALLY  ELASTIC  SCATTERING: 
THE  ELASTIC  CONTINUUM 


It  has  been  shown  in  Chapter  II  that  for  "well-behaved"  slowing-down 
operators  the  point  spectrum  of  the  wave  Boltzmann  operator  is  empty  and 
that  K  ^  C  belongs  to  the  resolvent  set.  Here  we  discuss  two  limiting 
processes  for  scattering  models  which  express  concretely  some  of  the 
ideas  presented  in  connection  with  almost-separable  mode  behavior,  and 
which  give  rise  to  continuous  spectra  for  k  ?  C. 

Duderstadt  [32,33]  discusses  the  effect  of  introducing  the 
macros copically  elastic  scattering  operator  of  Eq.  (2,25)  upon  the 
spectrum  of  the  wave  Boltzmann  operator.  Such  a  scattering  operator 
arises  in  the  theory  of  thermal  neutron  Bragg  scattering  in  poly- 
crystalline  materials.   In  Chapter  II  its  use  in  modeling  elastic 
scattering  from  heavy  nuclides  was  discussed  briefly. 

To  pursue  the  ideas  presented  there,  consider  for  simplicity  an 
isotropic  scattering  kernel  for  the  slowing-down  WBE  having  macro- 
scopically  elastic  and  nonelastic  terms 


Z^(E'->E)  =  SggCE')  a(E',E.AE)  +  E^^^(E'^E)  (C.l) 


where 


137 


138 


a(E',E,AE)  =  -^  , 


E'  -  AE  <  E  <  E- 


=  0 


otherwise 


(C.2) 


and  E   (E'->-E)  is  a  bounded  kernel.  Repeating  the  resolvent  set 
discussion  for  <  ^  C,   Eqs.  (2 .9)- (2 .15) ,  with  Eq.  (C.l)  as  the  kernel, 
we  find  that  the  complement  of  C  again  is  the  resolvent  set  for  any 
nonzero  AE,  since  both  kernels  are  bounded.  However  in  the  limit 
AE  =  0  Eq.  (C.l)  becomes  the  singular  kernel  of  Eq.  (2.25);  once  more 
repeating  the  resolvent  set  investigation  we  obtain  for  Eq.  (2.13) 


KE) 


=  f(E)iE^g(E)  (j)(E)  +  f(E)  j  j^sne^^'^^^  *(EOdE- 


g(E) 


(C.3) 


Defining 


AeCE) 


1  -  jfCE)  l3g(E) 


(C.4) 


Eq.  (C.3)  becomes 


<t.(E) 


_  f(E) 


A  (E) 


I  ^sne^^'^^^  ()>(E')dE'  +  g(E)  , 


K  g  C  . 


(C.5) 


For  values  of  k  (entering  through  f(E))  and  E  such  that  A  vanishes, 
Eq.  (C.5)  will  not  be  invertible  and  k   will  be  in  the  spectriom  of  the 
wave  Boltzmann  operator  [32,33].  For  a  particular  E,  Eq.  (C.4)  is  the 
dispersion  law  function  for  the  point  spectrum  of  the  monoenergetic 
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wave  transport  equation;  as  E  varies,  values  of  k  which  satisfy  the 
dispersion  law  will  sweep  out  a  line,  the  "elastic  continuum"  [32,33], 
in  the  complex  k   plane  exterior  to  the  transport  continuum  C.  For  a 
monoenergetic  source  this  additional  spectral  contribution  evidently 
will  give  rise  to  modes  which  appear  to  be  discrete  modes  having  the 
source  energy,  in  accordance  with  the  discussion  in  Chapter  II.  This 
appearance  of  additional  spectrum  only  in  the  limit  AE  =  0  is  not  sur- 
prising from  a  mathematical  point  of  view;  however,  physically  it  is 
evident  that  the  role  of  the  elastic  spectral  contribution  must  have 
been  assumed  for  AE  ;^  0  by  superimposed  continuum  eigenfunctions  for 
K  e  C.  Thus  we  must  be  cautious  in  regarding  continuum  modes  as  more 
attenuated  than  discrete  regular  modes,  since  superimposed  continuum 
modes  clearly  can  exhibit  highly  persistent  spatial  behavior. 

Another  simple  limiting  process  can  be  employed  to  model  the  same 
phenomenon,  making  use  of  the  elastic  continuum.  Define  the  kernel 

E  (E'^E)  =  BE  fE')  6(E'-E)  +  (1-6)  E  fE')  K(E'^E)  , 
s  s  s 

0  <_  6  <_  1  CC.6) 

so  that  B  acts  like  a  detector  energy  window,  determining  the  fraction 
of  neutrons  which  may  be  regarded  as  monoenergetic  after  one  collision. 
The  "elastic  continuum"  T  corresponding  to  this  model  is  shown 
schematically  in  Figure  C.l.  For  a  particular  energy  E  the  correspond- 
ing line  in  the  transport  continuum  is  denoted  by  < (E  ) ,  while  KqCEq) 
is  the  value  of  k  which  satisfies  the  elastic  dispersion  law  Ag(E^)  =  0. 
We  note  that  for  this  model  the  apparent  attenuation  and  wave  length  of 
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Figure  C.l    Schematic  diagram  of  the  "elastic  continuum" 
for  macroscopically  elastic  scattering. 
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the  elastic  mode  depend  on  the  value  of  6;  this  may  be  seen  by  using 
the  definition  of  f(E) 


£(E)  =  i  ln(2^)  (C.7) 


in  the  elastic  dispersion  function,  Eq.  (C.4)  so  that  Ag  =  0  may  be 
rearranged  to  give 


1  _   1  ^a+K, 

K    =    rr   I  ln( } 

2  se   ^a-K-^ 


=  1  8E^  ln(|il)  .  (C.8) 


For  <     not  too  near  the  continuum,  so  that  variations  in  the  logarithm 
e 

are  slow,  <  varies  almost  linearly  with  6,  as  shown  in  Figure  C.l. 
Evidently,  then,  the  observed  modal  propagation  constants  in  an  experi- 
ment would  depend  on  the  detector  energy  window. 
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